#################################################################################
### Title: Are People Willing to Trade Away Democracy for Desirable Outcomes? ###
### Authors: Jonathan A. Chu, Scott Williamson, Eddy S. F. Yeung              ###
### Content: Main analysis                                                    ###
### Date: September 20, 2025                                                  ###
#################################################################################

### Set-up ----
## Clean the working environment and set the working directory
rm(list = ls())
setwd("~/Desktop/democracy_tradeoff/replication") # change to your own working directory

## Load the required packages
library(tidyverse)
library(estimatr)
library(cregg)
library(expss)
library(cowplot)
library(grid)
library(gridExtra)
library(knitr)

## Read the cleaned and recoded datasets
df_US <- read.csv("df_US.csv")
df_JP <- read.csv("df_JP.csv")
df_EG <- read.csv("df_EG.csv")
df_IN <- read.csv("df_IN.csv")
df_IT <- read.csv("df_IT.csv")
df_TH <- read.csv("df_TH.csv")

## Merge the datasets
df_cj <- bind_rows(df_US, df_JP, df_EG, df_IN, df_IT, df_TH)

## Reorder the factors
# Leader selection
df_cj$`Leader Selection` <- 
  factor(df_cj$Leader.Selection,
         levels = c("Free and fair elections", "Unfair elections", "Unelected elites", "Hereditary succession", "Military coup"))

# Leader selection (binary)
df_cj$free_election <- ifelse(df_cj$Leader.Selection == "Free and fair elections", "Free and fair elections", "No free and fair elections")
df_cj$free_election <- 
  factor(df_cj$free_election,
         levels = c("Free and fair elections", "No free and fair elections"))

# Civil liberties
df_cj$`Civil Liberties` <- 
  factor(df_cj$Civil.Liberties,
         levels = c("Free", "Partially free", "Repressed"))

# Institutional checks
df_cj$`Institutional Checks` <- 
  factor(df_cj$Leader.Constraints,
         levels = c("Constrained", "Partially constrained", "Unconstrained"))

# National economy
df_cj$`National Economy` <- 
  factor(df_cj$National.Economy,
         levels = c("High income", "Middle income", "Low income"))

# Respondent wealth
df_cj$`Respondent Wealth` <- 
  factor(df_cj$Respondent.Wealth,
         levels = c("Wealthy", "Average", "Poor"))

# Public safety
df_cj$`Public Safety` <- 
  factor(df_cj$Public.Safety,
         levels = c("Very safe", "Somewhat safe", "Somewhat dangerous", "Very dangerous"))

# Corruption in politics
df_cj$`Corruption in Politics` <- 
  factor(df_cj$Corruption.in.Politics,
         levels = c("Rare", "Occasional", "Prevalent"))

# Health care
df_cj$`Health Care` <- 
  factor(df_cj$Health.Care,
         levels = c("Mostly accessible", "For the privileged"))

# Minority treatment
df_cj$`Minority Treatment` <- 
  factor(df_cj$Minority.Treatment,
         levels = c("Fairly treated", "Sometimes unfair", "Mostly unfair"))

# Respondent identity
df_cj$`Respondent Identity` <- 
  factor(df_cj$Respondent.Identity,
         levels = c("Majority", "Second largest", "Minority"))

# Country
df_cj$country <- 
  factor(df_cj$country, 
         levels = c("US", "JP", "EG", "IN", "IT", "TH"),
         labels = c("United States", "Japan", "Egypt", "India", "Italy", "Thailand"))

### Plot the AMCEs in aggregate sample (Figure 1) ----
## Bold feature labels in plots
bph <- c('bold', rep('plain', 5), 'bold', rep('plain', 3),
         'bold', rep('plain', 3), 'bold', rep('plain', 3), 
         'bold', rep('plain', 3), 'bold', rep('plain', 4),
         'bold', rep('plain', 3), 'bold', rep('plain', 2), 
         'bold', rep('plain', 3), 'bold', rep('plain', 3)) %>%
  rev() # reverse coding

## Create a function that indicates which estimates are statistically significant
## at the 5% level after using the BH procedure to adjust for multiple comparisons
sig1_fun <- function(data){
  data <- data %>% 
    mutate(sig = case_when(
      (p.adjust(p, method = "BH") < 0.05) == TRUE ~ 1, 
      (p.adjust(p, method = "BH") < 0.05) == FALSE ~ 0)) %>% 
    mutate(sig = factor(sig, levels = c(0, 1)))
  return(data)
}

## Estimate the AMCEs
amces <- cj(data = df_cj,
            formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
              `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`,
            level_order = "descending",
            id = ~id,
            estimate = "amce")
amces <- sig1_fun(amces)
amces$sig[is.na(amces$sig)] <- 0

## Visualize the AMCEs
p <- plot(amces, legend_pos = "none", size = 1.5) + 
  aes(color = sig) + 
  xlab("Effect on Pr(preferring to be born and grow up in the country)")
p <- p + scale_color_manual(values = c("grey60", "black"))
p <- p + theme_bw(base_size = 10, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 11, colour = "black",  hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 10, colour = "black", hjust = 1, vjust = .5, face = bph),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 11),
    axis.title.y = element_text(size = 11, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 10),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()
  )
p_amce <- p + coord_cartesian(xlim = c(-0.225, 0.025))
p_amce
ggsave("Figure 1.pdf", width = 8, height = 6)

## Present results in table form (Table S5)
amces <- subset(amces, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"))
amces$N <- 36900
kable(amces, digits = 3, "latex")

### Plot the marginal means by country (Figures S1-S6) ----
## Egypt (Figure S1)
mm <- cj(data = subset(df_cj, country == "Egypt"),
         formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
           `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
           `Health Care` + `Minority Treatment` + `Respondent Identity`,
         level_order = "descending",
         id = ~id,
         estimate = "mm")
p <- plot(mm, legend_pos = "none", size = 1.5) + 
  aes(color = as.factor(1)) +
  xlab("Pr(preferring to be born and grow up in the country)")
p <- p + scale_color_manual(values = c("black"))
p <- p + theme_bw(base_size = 10, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 11, colour = "black",  hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 10, colour = "black", hjust = 1, vjust = .5, face = bph),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 11),
    axis.title.y = element_text(size = 11, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 10),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()
  )
p_mm <- p + coord_cartesian(xlim = c(0.3, 0.7)) +
  geom_vline(xintercept = 0.5, color = "gray75", linetype = "dashed") +
  theme(plot.title = element_text(hjust = 0.5))
p_mm
ggsave("Figure S1.pdf", width = 8, height = 6)

# Present results in table form (Table S18)
mm_table <- subset(mm, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"))
mm_table$N <- 6048
kable(mm_table, digits = 3, "latex")

## India (Figure S2)
mm <- cj(data = subset(df_cj, country == "India"),
         formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
           `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
           `Health Care` + `Minority Treatment` + `Respondent Identity`,
         level_order = "descending",
         id = ~id,
         estimate = "mm")
p <- plot(mm, legend_pos = "none", size = 1.5) + 
  aes(color = as.factor(1)) +
  xlab("Pr(preferring to be born and grow up in the country)")
p <- p + scale_color_manual(values = c("black"))
p <- p + theme_bw(base_size = 10, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 11, colour = "black",  hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 10, colour = "black", hjust = 1, vjust = .5, face = bph),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 11),
    axis.title.y = element_text(size = 11, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 10),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()
  )
p_mm <- p + coord_cartesian(xlim = c(0.3, 0.7)) +
  geom_vline(xintercept = 0.5, color = "gray75", linetype = "dashed") +
  theme(plot.title = element_text(hjust = 0.5))
p_mm
ggsave("Figure S2.pdf", width = 8, height = 6)

# Present results in table form (Table S19)
mm_table <- subset(mm, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"))
mm_table$N <- 6132
kable(mm_table, digits = 3, "latex")

## Italy (Figure S3)
mm <- cj(data = subset(df_cj, country == "Italy"),
         formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
           `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
           `Health Care` + `Minority Treatment` + `Respondent Identity`,
         level_order = "descending",
         id = ~id,
         estimate = "mm")
p <- plot(mm, legend_pos = "none", size = 1.5) + 
  aes(color = as.factor(1)) +
  xlab("Pr(preferring to be born and grow up in the country)")
p <- p + scale_color_manual(values = c("black"))
p <- p + theme_bw(base_size = 10, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 11, colour = "black",  hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 10, colour = "black", hjust = 1, vjust = .5, face = bph),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 11),
    axis.title.y = element_text(size = 11, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 10),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()
  )
p_mm <- p + coord_cartesian(xlim = c(0.3, 0.7)) +
  geom_vline(xintercept = 0.5, color = "gray75", linetype = "dashed") +
  theme(plot.title = element_text(hjust = 0.5))
p_mm
ggsave("Figure S3.pdf", width = 8, height = 6)

# Present results in table form (Table S20)
mm_table <- subset(mm, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"))
mm_table$N <- 6282
kable(mm_table, digits = 3, "latex")

## Japan (Figure S4)
mm <- cj(data = subset(df_cj, country == "Japan"),
         formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
           `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
           `Health Care` + `Minority Treatment` + `Respondent Identity`,
         level_order = "descending",
         id = ~id,
         estimate = "mm")
p <- plot(mm, legend_pos = "none", size = 1.5) + 
  aes(color = as.factor(1)) +
  xlab("Pr(preferring to be born and grow up in the country)")
p <- p + scale_color_manual(values = c("black"))
p <- p + theme_bw(base_size = 10, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 11, colour = "black",  hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 10, colour = "black", hjust = 1, vjust = .5, face = bph),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 11),
    axis.title.y = element_text(size = 11, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 10),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()
  )
p_mm <- p + coord_cartesian(xlim = c(0.3, 0.7)) +
  geom_vline(xintercept = 0.5, color = "gray75", linetype = "dashed") +
  theme(plot.title = element_text(hjust = 0.5))
p_mm
ggsave("Figure S4.pdf", width = 8, height = 6)

# Present results in table form (Table S21)
mm_table <- subset(mm, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"))
mm_table$N <- 6072
kable(mm_table, digits = 3, "latex")

## Thailand (Figure S5)
mm <- cj(data = subset(df_cj, country == "Thailand"),
         formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
           `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
           `Health Care` + `Minority Treatment` + `Respondent Identity`,
         level_order = "descending",
         id = ~id,
         estimate = "mm")
p <- plot(mm, legend_pos = "none", size = 1.5) + 
  aes(color = as.factor(1)) +
  xlab("Pr(preferring to be born and grow up in the country)")
p <- p + scale_color_manual(values = c("black"))
p <- p + theme_bw(base_size = 10, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 11, colour = "black",  hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 10, colour = "black", hjust = 1, vjust = .5, face = bph),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 11),
    axis.title.y = element_text(size = 11, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 10),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()
  )
p_mm <- p + coord_cartesian(xlim = c(0.3, 0.7)) +
  geom_vline(xintercept = 0.5, color = "gray75", linetype = "dashed") +
  theme(plot.title = element_text(hjust = 0.5))
p_mm
ggsave("Figure S5.pdf", width = 8, height = 6)

# Present results in table form (Table S22)
mm_table <- subset(mm, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"))
mm_table$N <- 6222
kable(mm_table, digits = 3, "latex")

## United States (Figure S6)
mm <- cj(data = subset(df_cj, country == "United States"),
         formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
           `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
           `Health Care` + `Minority Treatment` + `Respondent Identity`,
         level_order = "descending",
         id = ~id,
         estimate = "mm")
p <- plot(mm, legend_pos = "none", size = 1.5) + 
  aes(color = as.factor(1)) +
  xlab("Pr(preferring to be born and grow up in the country)")
p <- p + scale_color_manual(values = c("black"))
p <- p + theme_bw(base_size = 10, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 11, colour = "black",  hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 10, colour = "black", hjust = 1, vjust = .5, face = bph),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 11),
    axis.title.y = element_text(size = 11, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 10),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()
  )
p_mm <- p + coord_cartesian(xlim = c(0.3, 0.7)) +
  geom_vline(xintercept = 0.5, color = "gray75", linetype = "dashed") +
  theme(plot.title = element_text(hjust = 0.5))
p_mm
ggsave("Figure S6.pdf", width = 8, height = 6)

# Present results in table form (Table S23)
mm_table <- subset(mm, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"))
mm_table$N <- 6144
kable(mm_table, digits = 3, "latex")

### Plot the AMCEs for individual countries (Figure 2) ----
## Estimate the AMCEs
amces <- 
  cj(data = df_cj, 
     formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
       `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
       `Health Care` + `Minority Treatment` + `Respondent Identity`,
     level_order = "descending",
     id = ~id,
     estimate = "amce",
     by = ~country)
amces <- sig1_fun(amces)

## Visualize the AMCEs
amces$estimate[amces$BY == "Japan" & amces$level == "Somewhat dangerous"] <- NA
p <- plot(amces, legend_pos = "none", size = 1) + 
  ggplot2::facet_wrap(~BY, ncol = 3) +
  aes(color = sig) + 
  xlab("Effect on Pr(preferring to be born and grow up in the country)")
p <- p + scale_color_manual(values = c("grey60", "black"))
p <- p + theme_bw(base_size = 8, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 8, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 8, colour = "black", hjust = 1, vjust = .5, face = bph),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 8)
  )
p <- p + coord_cartesian(xlim = c(-0.3, 0.05))
p
ggsave("Figure 2.pdf", width = 6.5, height = 8.5)

## Present results in table form (Tables S6-S11)
amces_EG <- subset(amces, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"),
                   BY == "Egypt")
amces_EG$N <- nrow(df_EG)
kable(amces_EG, digits = 3, "latex", row.names = F)

amces_IN <- subset(amces, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"),
                   BY == "India")
amces_IN$N <- nrow(df_IN)
kable(amces_IN, digits = 3, "latex", row.names = F)

amces_IT <- subset(amces, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"),
                   BY == "Italy")
amces_IT$N <- nrow(df_IT)
kable(amces_IT, digits = 3, "latex", row.names = F)

amces_JP <- subset(amces, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"),
                   BY == "Japan")
amces_JP$N <- nrow(df_JP)
kable(amces_JP, digits = 3, "latex", row.names = F)

amces_TH <- subset(amces, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"),
                   BY == "Thailand")
amces_TH$N <- nrow(df_TH)
kable(amces_TH, digits = 3, "latex", row.names = F)

amces_US <- subset(amces, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"),
                   BY == "United States")
amces_US$N <- nrow(df_US)
kable(amces_US, digits = 3, "latex", row.names = F)

### Plot differences in marginal means by subgroup (Figure S10) ----
## By age
df_cj$age_bin <- factor(df_cj$age_bin, levels = c("Older", "Younger Than 40"))
diff_mms <- cj(data = df_cj, 
               formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
                 `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
                 `Health Care` + `Minority Treatment` + `Respondent Identity`,
               level_order = "descending",
               id = ~id,
               estimate = "mm_diff",
               by = ~age_bin)
diff_mms_age <- sig1_fun(diff_mms)

## By gender
df_cj$gender_bin <- factor(df_cj$gender_bin, levels = c("Male", "Female"))
diff_mms <- cj(data = df_cj, 
               formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
                 `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
                 `Health Care` + `Minority Treatment` + `Respondent Identity`,
               level_order = "descending",
               id = ~id,
               estimate = "mm_diff",
               by = ~gender_bin)
diff_mms_gender <- sig1_fun(diff_mms)

## By education
df_cj$edu_bin <- factor(df_cj$edu_bin, levels = c("No College", "College"))
diff_mms <- cj(data = df_cj, 
               formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
                 `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
                 `Health Care` + `Minority Treatment` + `Respondent Identity`,
               level_order = "descending",
               id = ~id,
               estimate = "mm_diff",
               by = ~edu_bin)
diff_mms_edu <- sig1_fun(diff_mms)

## By socioeconomic status
df_cj$SES <- factor(df_cj$SES, levels = c("Low SES", "High SES"))
diff_mms <- cj(data = df_cj, 
               formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
                 `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
                 `Health Care` + `Minority Treatment` + `Respondent Identity`,
               level_order = "descending",
               id = ~id,
               estimate = "mm_diff",
               by = ~SES)
diff_mms_SES <- sig1_fun(diff_mms)

## By minority status
df_cj$minority_bin <- factor(df_cj$minority_bin, levels = c("Non-Minority", "Minority"))
diff_mms <- cj(data = df_cj, 
               formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
                 `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
                 `Health Care` + `Minority Treatment` + `Respondent Identity`,
               level_order = "descending",
               id = ~id,
               estimate = "mm_diff",
               by = ~minority_bin)
diff_mms_minority <- sig1_fun(diff_mms)

## By self-reported political ideology
df_cj$ideo_bin <- factor(df_cj$ideo_bin, levels = c("Leftwing", "Rightwing"))
diff_mms <- cj(data = df_cj, 
               formula = selected ~ `Leader Selection` + `Civil Liberties` + `Institutional Checks` +
                 `National Economy` + `Respondent Wealth` + `Public Safety` + `Corruption in Politics` + 
                 `Health Care` + `Minority Treatment` + `Respondent Identity`,
               level_order = "descending",
               id = ~id,
               estimate = "mm_diff",
               by = ~ideo_bin)
diff_mms_ideo <- sig1_fun(diff_mms)

## Combine the difference-in-MMs plots
combined <- bind_rows(diff_mms_age, diff_mms_gender, diff_mms_edu,
                      diff_mms_SES, diff_mms_minority, diff_mms_ideo)
combined$BY <- 
  factor(combined$BY,
         levels = c("Younger Than 40 - Older", "Female - Male", "College - No College",
                    "High SES - Low SES", "Minority - Non-Minority", "Rightwing - Leftwing"))
p <- plot(combined, legend_pos = "none", size = 1) + 
  ggplot2::facet_wrap(~BY, ncol = 3) +
  aes(color = sig) +
  xlab("Difference in marginal means")
p <- p + scale_color_manual(values = c("grey60", "black"))
p <- p + theme_bw(base_size = 8, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 8, colour = "black",  hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 8, colour = "black", hjust = 1, vjust = .5, face = bph),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 8)
  )
p <- p + coord_cartesian(xlim = c(-0.07, 0.07))
p
ggsave("Figure S10.pdf", width = 6.5, height = 8.5)

## Present results in table form (Tables S26-S31)
diff_mms_age <- subset(diff_mms_age, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"))
table(df_cj$age_bin)
diff_mms_age$N <- 36120
kable(diff_mms_age, digits = 3, "latex")

diff_mms_gender <- subset(diff_mms_gender, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"))
table(df_cj$gender_bin)
diff_mms_gender$N <- 36510
kable(diff_mms_gender, digits = 3, "latex")

diff_mms_edu <- subset(diff_mms_edu, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"))
table(df_cj$edu_bin)
diff_mms_edu$N <- 36900
kable(diff_mms_edu, digits = 3, "latex")

diff_mms_SES <- subset(diff_mms_SES, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"))
table(df_cj$SES)
diff_mms_SES$N <- 36894
kable(diff_mms_SES, digits = 3, "latex")

diff_mms_minority <- subset(diff_mms_minority, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"))
table(df_cj$minority_bin)
diff_mms_minority$N <- 36900
kable(diff_mms_minority, digits = 3, "latex")

diff_mms_ideo <- subset(diff_mms_ideo, select = c("feature", "level", "estimate", "std.error", "p", "lower", "upper"))
table(df_cj$ideo_bin)
diff_mms_ideo$N <- 18486
kable(diff_mms_ideo, digits = 3, "latex")

### Analyze interaction effects between leader selection and public safety (Figure 3A) ----
## Create a new interaction variable
df_cj$election_safety <- interaction(df_cj$free_election, df_cj$`Public Safety`, sep = " × ")
df_cj$election_safety <- relevel(df_cj$election_safety, ref = "Free and fair elections × Very dangerous")

## Estimate the average component interaction effects
acie_election_safety_est <- 
  lm_robust(selected ~ election_safety + `Civil Liberties` + `Institutional Checks` +
              `Corruption in Politics` + `National Economy` + `Respondent Wealth` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, clusters = id, fixed_effects = country)
acie_election_safety <- acie_election_safety_est %>% 
  tidy() %>% 
  add_row(term = "Free and fair elections × Very dangerous", estimate = 0) %>% 
  slice(1:7, 23)
acie_election_safety$term <- 
  factor(acie_election_safety$term,
         levels = c("election_safetyNo free and fair elections × Very dangerous",
                    "election_safetyNo free and fair elections × Somewhat dangerous",
                    "election_safetyNo free and fair elections × Somewhat safe",
                    "election_safetyNo free and fair elections × Very safe",
                    "Free and fair elections × Very dangerous",
                    "election_safetyFree and fair elections × Somewhat dangerous",
                    "election_safetyFree and fair elections × Somewhat safe",
                    "election_safetyFree and fair elections × Very safe"),
         labels = c("No free and fair elections ×\nVery dangerous",
                    "No free and fair elections ×\nSomewhat dangerous",
                    "No free and fair elections ×\nSomewhat safe",
                    "No free and fair elections ×\nVery safe",
                    "Free and fair elections ×\nVery dangerous",
                    "Free and fair elections ×\nSomewhat dangerous",
                    "Free and fair elections ×\nSomewhat safe",
                    "Free and fair elections ×\nVery safe"))
acie_election_safety$color <- ifelse(acie_election_safety$estimate == 0, 0, 1)

p1 <- ggplot(acie_election_safety, aes(x = estimate, y = term, color = factor(color))) +
  geom_point(size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0) +
  xlab("Effect on Pr(preferring to be born\nand grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Free and Fair Elections vs. Public Safety") +
  coord_cartesian(xlim = c(-.25, .25))
p1 <- p1 + scale_color_manual(values = c("grey60", "black"))
p1 <- p1 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 11, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white")
  )

### Analyze interaction effects between leader selection and respondent wealth (Figure 3B) ----
## Create a new interaction variable
df_cj$election_wealth <- interaction(df_cj$free_election, df_cj$`Respondent Wealth`, sep = " × ")
df_cj$election_wealth <- relevel(df_cj$election_wealth, ref = "Free and fair elections × Poor")

## Estimate the average component interaction effects
acie_election_wealth_est <- 
  lm_robust(selected ~ election_wealth + `Civil Liberties` + `Institutional Checks` +
              `Corruption in Politics` + `National Economy` + `Public Safety` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, clusters = id, fixed_effects = country)
acie_election_wealth <- acie_election_wealth_est %>% 
  tidy() %>% 
  add_row(term = "Free and fair elections × Poor", estimate = 0) %>% 
  slice(1:5, 22)
acie_election_wealth$term <- 
  factor(acie_election_wealth$term,
         levels = c("election_wealthNo free and fair elections × Poor",
                    "election_wealthNo free and fair elections × Average",
                    "election_wealthNo free and fair elections × Wealthy",
                    "Free and fair elections × Poor",
                    "election_wealthFree and fair elections × Average",
                    "election_wealthFree and fair elections × Wealthy"),
         labels = c("No free and fair elections ×\nPoor",
                    "No free and fair elections ×\nAverage",
                    "No free and fair elections ×\nWealthy",
                    "Free and fair elections ×\nPoor",
                    "Free and fair elections ×\nAverage",
                    "Free and fair elections ×\nWealthy"))
acie_election_wealth$color <- ifelse(acie_election_wealth$estimate == 0, 0, 1)

p2 <- ggplot(acie_election_wealth, aes(x = estimate, y = term, color = factor(color))) +
  geom_point(size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0) +
  xlab("Effect on Pr(preferring to be born\nand grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Free and Fair Elections vs. Respondent Wealth") +
  coord_cartesian(xlim = c(-.25, .25))
p2 <- p2 + scale_color_manual(values = c("grey60", "black"))
p2 <- p2 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 11, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white")
  )

plot_grid(p1, p2, 
          labels = "AUTO", ncol = 2, 
          label_fontfamily = "Times")
ggsave("Figure 3.pdf", width = 10, height = 6)

## Present results in table form (Tables S12-S13)
acies_safety <- subset(acie_election_safety, select = c("term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
acies_safety$N <- 36900
kable(acies_safety, digits = 3, "latex")

acies_wealth <- subset(acie_election_wealth, select = c("term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
acies_wealth$N <- 36900
kable(acies_wealth, digits = 3, "latex")

### Analyze interaction effects between leader selection and national economy (Figure S11) ----
## Create a new interaction variable
df_cj$election_econ <- interaction(df_cj$free_election, df_cj$`National Economy`, sep = " × ")
df_cj$election_econ <- 
  relevel(df_cj$election_econ, ref = "Free and fair elections × Low income")

## Estimate the average component interaction effects
acie_election_econ_est <- 
  lm_robust(selected ~ election_econ + `Civil Liberties` + `Institutional Checks` +
              `Corruption in Politics` + `Respondent Wealth` + `Public Safety` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, clusters = id, fixed_effects = country)
acie_election_econ <- acie_election_econ_est %>% 
  tidy() %>% 
  add_row(term = "Free and fair elections × Low income", estimate = 0) %>% 
  slice(1:5, 22)
acie_election_econ$term <- 
  factor(acie_election_econ$term,
         levels = c("election_econNo free and fair elections × Low income",
                    "election_econNo free and fair elections × Middle income",
                    "election_econNo free and fair elections × High income",
                    "Free and fair elections × Low income",
                    "election_econFree and fair elections × Middle income",
                    "election_econFree and fair elections × High income"),
         labels = c("No free and fair elections ×\nLow income",
                    "No free and fair elections ×\nMiddle income",
                    "No free and fair elections ×\nHigh income",
                    "Free and fair elections ×\nLow income",
                    "Free and fair elections ×\nMiddle income",
                    "Free and fair elections ×\nHigh income"))
acie_election_econ$color <- ifelse(acie_election_econ$estimate == 0, 0, 1)

p3 <- ggplot(acie_election_econ, aes(x = estimate, y = term, color = factor(color))) +
  geom_point(size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0) +
  xlab("Effect on Pr(preferring to be born\nand grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Free and Fair Elections vs. National Economy") +
  coord_cartesian(xlim = c(-.25, .25))
p3 <- p3 + scale_color_manual(values = c("grey60", "black"))
p3 <- p3 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 11, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white")
  )
p3
ggsave("Figure S11.pdf", width = 6, height = 6)

## Present results in table form (Table S32)
acies_econ <- subset(acie_election_econ, select = c("term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
acies_econ$N <- 36900
kable(acies_econ, digits = 3, "latex")

### Analyze interaction effects between leader selection and health care (Figure S14) ----
## Create a new interaction variable
df_cj$election_heatlh <- interaction(df_cj$free_election, df_cj$`Health Care`, sep = " × ")
df_cj$election_heatlh <- 
  relevel(df_cj$election_heatlh, ref = "Free and fair elections × For the privileged")

## Estimate the average component interaction effects
acie_election_health_est <- 
  lm_robust(selected ~ election_heatlh + `Civil Liberties` + `Institutional Checks` +
              `Corruption in Politics` + `National Economy` + `Respondent Wealth` + 
              `Public Safety` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, clusters = id, fixed_effects = country)
acie_election_health <- acie_election_health_est %>% 
  tidy() %>% 
  add_row(term = "Free and fair elections × For the privileged", estimate = 0) %>% 
  slice(1:3, 21)
acie_election_health$term <- 
  factor(acie_election_health$term,
         levels = c("election_heatlhNo free and fair elections × For the privileged",
                    "election_heatlhNo free and fair elections × Mostly accessible",
                    "Free and fair elections × For the privileged",
                    "election_heatlhFree and fair elections × Mostly accessible"),
         labels = c("No free and fair elections ×\nHealth care for the privileged",
                    "No free and fair elections ×\nHealth care mostly accessible",
                    "Free and fair elections ×\nHealth care for the privileged",
                    "Free and fair elections ×\nHealth care mostly accessible"))
acie_election_health$color <- ifelse(acie_election_health$estimate == 0, 0, 1)

p4 <- ggplot(acie_election_health, aes(x = estimate, y = term, color = factor(color))) +
  geom_point(size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0) +
  xlab("Effect on Pr(preferring to be born\nand grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Free and Fair Elections vs. Health Care") +
  coord_cartesian(xlim = c(-.25, .25))
p4 <- p4 + scale_color_manual(values = c("grey60", "black"))
p4 <- p4 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 11, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white")
  )
p4
ggsave("Figure S14.pdf", width = 6, height = 6)

## Present results in table form (Table S35)
acies_health <- subset(acie_election_health, select = c("term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
acies_health$N <- 36900
kable(acies_health, digits = 3, "latex")

### Analyze interaction effects between civil liberties and public safety (Figure 4A) ----
## Create a new interaction variable
df_cj$civil_safety <- interaction(df_cj$`Civil Liberties`, df_cj$`Public Safety`, sep = " × ")
df_cj$civil_safety <- relevel(df_cj$civil_safety, ref = "Free × Very dangerous")

## Estimate the average component interaction effects
acie_civil_safety_est <- 
  lm_robust(selected ~ civil_safety + `Leader Selection` + `Institutional Checks` +
              `Corruption in Politics` + `National Economy` + `Respondent Wealth` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, clusters = id, fixed_effects = country)
acie_civil_safety <- acie_civil_safety_est %>% 
  tidy() %>% 
  add_row(term = "Free × Very dangerous", estimate = 0) %>% 
  slice(1:11, 29)
acie_civil_safety$term <- 
  factor(acie_civil_safety$term,
         levels = c("civil_safetyRepressed × Very dangerous",
                    "civil_safetyRepressed × Somewhat dangerous",
                    "civil_safetyRepressed × Somewhat safe",
                    "civil_safetyRepressed × Very safe",
                    "civil_safetyPartially free × Very dangerous",
                    "civil_safetyPartially free × Somewhat dangerous",
                    "civil_safetyPartially free × Somewhat safe",
                    "civil_safetyPartially free × Very safe",
                    "Free × Very dangerous",
                    "civil_safetyFree × Somewhat dangerous",
                    "civil_safetyFree × Somewhat safe",
                    "civil_safetyFree × Very safe"),
         labels = c("Repressed ×\nVery dangerous",
                    "Repressed ×\nSomewhat dangerous",
                    "Repressed ×\nSomewhat safe",
                    "Repressed ×\nVery safe",
                    "Partially free ×\nVery dangerous",
                    "Partially free ×\nSomewhat dangerous",
                    "Partially free ×\nSomewhat safe",
                    "Partially free ×\nVery safe",
                    "Free ×\nVery dangerous",
                    "Free ×\nSomewhat dangerous",
                    "Free ×\nSomewhat safe",
                    "Free ×\nVery safe"))
acie_civil_safety$color <- ifelse(acie_civil_safety$estimate == 0, 0, 1)

p1 <- ggplot(acie_civil_safety, aes(x = estimate, y = term, color = factor(color))) +
  geom_point(size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0) +
  xlab("Effect on Pr(preferring to be born\nand grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Civil Liberties vs. Public Safety") +
  coord_cartesian(xlim = c(-.25, .25))
p1 <- p1 + scale_color_manual(values = c("grey60", "black"))
p1 <- p1 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 11, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white")
  )

### Analyze interaction effects between civil liberties and respondent wealth (Figure 4B) ----
## Create a new interaction variable
df_cj$civil_wealth <- interaction(df_cj$`Civil Liberties`, df_cj$`Respondent Wealth`, sep = " × ")
df_cj$civil_wealth <- relevel(df_cj$civil_wealth, ref = "Free × Poor")

## Estimate the average component interaction effects
acie_civil_wealth_est <- 
  lm_robust(selected ~ civil_wealth + `Leader Selection` + `Institutional Checks` +
              `Corruption in Politics` + `National Economy` + `Public Safety` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, clusters = id, fixed_effects = country)
acie_civil_wealth <- acie_civil_wealth_est %>% 
  tidy() %>% 
  add_row(term = "Free × Poor", estimate = 0) %>% 
  slice(1:8, 27)
acie_civil_wealth$term <- 
  factor(acie_civil_wealth$term,
         levels = c("civil_wealthRepressed × Poor",
                    "civil_wealthRepressed × Average",
                    "civil_wealthRepressed × Wealthy",
                    "civil_wealthPartially free × Poor",
                    "civil_wealthPartially free × Average",
                    "civil_wealthPartially free × Wealthy",
                    "Free × Poor",
                    "civil_wealthFree × Average",
                    "civil_wealthFree × Wealthy"),
         labels = c("Repressed ×\nPoor",
                    "Repressed ×\nAverage",
                    "Repressed ×\nWealthy",
                    "Partially free ×\nPoor",
                    "Partially free ×\nAverage",
                    "Partially free ×\nWealthy",
                    "Free ×\nPoor",
                    "Free ×\nAverage",
                    "Free ×\nWealthy"))
acie_civil_wealth$color <- ifelse(acie_civil_wealth$estimate == 0, 0, 1)

p2 <- ggplot(acie_civil_wealth, aes(x = estimate, y = term, color = factor(color))) +
  geom_point(size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0) +
  xlab("Effect on Pr(preferring to be born\nand grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Civil Liberties vs. Respondent Wealth") +
  coord_cartesian(xlim = c(-.25, .25))
p2 <- p2 + scale_color_manual(values = c("grey60", "black"))
p2 <- p2 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 11, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white")
  )

plot_grid(p1, p2, 
          labels = "AUTO", ncol = 2, 
          label_fontfamily = "Times")
ggsave("Figure 4.pdf", width = 10, height = 6)

## Present results in table form (Tables S14-S15)
acies_safety <- subset(acie_civil_safety, select = c("term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
acies_safety$N <- 36900
kable(acies_safety, digits = 3, "latex")

acies_wealth <- subset(acie_civil_wealth, select = c("term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
acies_wealth$N <- 36900
kable(acies_wealth, digits = 3, "latex")

### Analyze interaction effects between civil liberties and national economy (Figure S12) ----
## Create a new interaction variable
df_cj$civil_econ <- interaction(df_cj$`Civil Liberties`, df_cj$`National Economy`, sep = " × ")
df_cj$civil_econ <- relevel(df_cj$civil_econ, ref = "Free × Low income")

## Estimate the average component interaction effects
acie_civil_econ_est <- 
  lm_robust(selected ~ civil_econ + `Leader Selection` + `Institutional Checks` +
              `Corruption in Politics` + `Respondent Wealth` + `Public Safety` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, clusters = id, fixed_effects = country)
acie_civil_econ <- acie_civil_econ_est %>% 
  tidy() %>% 
  add_row(term = "Free × Low income", estimate = 0) %>% 
  slice(1:8, 27)
acie_civil_econ$term <- 
  factor(acie_civil_econ$term,
         levels = c("civil_econRepressed × Low income",
                    "civil_econRepressed × Middle income",
                    "civil_econRepressed × High income",
                    "civil_econPartially free × Low income",
                    "civil_econPartially free × Middle income",
                    "civil_econPartially free × High income",
                    "Free × Low income",
                    "civil_econFree × Middle income",
                    "civil_econFree × High income"),
         labels = c("Repressed ×\nLow income",
                    "Repressed ×\nMiddle income",
                    "Repressed ×\nHigh income",
                    "Partially free ×\nLow income",
                    "Partially free ×\nMiddle income",
                    "Partially free ×\nHigh income",
                    "Free ×\nLow income",
                    "Free ×\nMiddle income",
                    "Free ×\nHigh income"))
acie_civil_econ$color <- ifelse(acie_civil_econ$estimate == 0, 0, 1)

p3 <- ggplot(acie_civil_econ, aes(x = estimate, y = term, color = factor(color))) +
  geom_point(size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0) +
  xlab("Effect on Pr(preferring to be born\nand grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Civil Liberties vs. National Economy") +
  coord_cartesian(xlim = c(-.25, .25))
p3 <- p3 + scale_color_manual(values = c("grey60", "black"))
p3 <- p3 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 11, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white")
  )
p3
ggsave("Figure S12.pdf", width = 6, height = 6)

## Present results in table form (Table S33)
acies_econ <- subset(acie_civil_econ, select = c("term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
acies_econ$N <- 36900
kable(acies_econ, digits = 3, "latex")

### Analyze interaction effects between civil liberties and health care (Figure S15) ----
## Create a new interaction variable
df_cj$civil_health <- interaction(df_cj$`Civil Liberties`, df_cj$`Health Care`, sep = " × ")
df_cj$civil_health <- relevel(df_cj$civil_health, ref = "Free × For the privileged")

## Estimate the average component interaction effects
acie_civil_health_est <- 
  lm_robust(selected ~ civil_health + `Leader Selection` + `Institutional Checks` +
              `Corruption in Politics` + `National Economy` + `Respondent Wealth` + 
              `Public Safety` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, clusters = id, fixed_effects = country)
acie_civil_health <- acie_civil_health_est %>% 
  tidy() %>% 
  add_row(term = "Free × For the privileged", estimate = 0) %>% 
  slice(1:5, 25)
acie_civil_health$term <- 
  factor(acie_civil_health$term,
         levels = c("civil_healthRepressed × For the privileged",
                    "civil_healthRepressed × Mostly accessible",
                    "civil_healthPartially free × For the privileged",
                    "civil_healthPartially free × Mostly accessible",
                    "Free × For the privileged",
                    "civil_healthFree × Mostly accessible"),
         labels = c("Repressed ×\nHealth care for the privileged",
                    "Repressed ×\nHealth care mostly accessible",
                    "Partially free ×\nHealth care for the privileged",
                    "Partially free ×\nHealth care mostly accessible",
                    "Free ×\nHealth care for the privileged",
                    "Free ×\nHealth care mostly accessible"))
acie_civil_health$color <- ifelse(acie_civil_health$estimate == 0, 0, 1)

p4 <- ggplot(acie_civil_health, aes(x = estimate, y = term, color = factor(color))) +
  geom_point(size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0) +
  xlab("Effect on Pr(preferring to be born\nand grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Civil Liberties vs. Health Care") +
  coord_cartesian(xlim = c(-.25, .25))
p4 <- p4 + scale_color_manual(values = c("grey60", "black"))
p4 <- p4 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 11, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white")
  )
p4
ggsave("Figure S15.pdf", width = 6, height = 6)

## Present results in table form (Table S36)
acies_health <- subset(acie_civil_health, select = c("term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
acies_health$N <- 36900
kable(acies_health, digits = 3, "latex")

### Analyze interaction effects between institutional checks and public safety (Figure 5A) ----
## Create a new interaction variable
df_cj$constr_safety <- interaction(df_cj$`Institutional Checks`, df_cj$`Public Safety`, sep = " × ")
df_cj$constr_safety <- relevel(df_cj$constr_safety, ref = "Constrained × Very dangerous")

## Estimate the average component interaction effects
acie_constr_safety_est <- 
  lm_robust(selected ~ constr_safety + `Leader Selection` + `Civil Liberties` +
              `Corruption in Politics` + `National Economy` + `Respondent Wealth` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, clusters = id, fixed_effects = country)
acie_constr_safety <- acie_constr_safety_est %>% 
  tidy() %>% 
  add_row(term = "Constrained × Very dangerous", estimate = 0) %>% 
  slice(1:11, 29)
acie_constr_safety$term <- 
  factor(acie_constr_safety$term,
         levels = c("constr_safetyUnconstrained × Very dangerous",
                    "constr_safetyUnconstrained × Somewhat dangerous",
                    "constr_safetyUnconstrained × Somewhat safe",
                    "constr_safetyUnconstrained × Very safe",
                    "constr_safetyPartially constrained × Very dangerous",
                    "constr_safetyPartially constrained × Somewhat dangerous",
                    "constr_safetyPartially constrained × Somewhat safe",
                    "constr_safetyPartially constrained × Very safe",
                    "Constrained × Very dangerous",
                    "constr_safetyConstrained × Somewhat dangerous",
                    "constr_safetyConstrained × Somewhat safe",
                    "constr_safetyConstrained × Very safe"),
         labels = c("Unconstrained ×\nVery dangerous",
                    "Unconstrained ×\nSomewhat dangerous",
                    "Unconstrained ×\nSomewhat safe",
                    "Unconstrained ×\nVery safe",
                    "Partially constrained ×\nVery dangerous",
                    "Partially constrained ×\nSomewhat dangerous",
                    "Partially constrained ×\nSomewhat safe",
                    "Partially constrained ×\nVery safe",
                    "Constrained ×\nVery dangerous",
                    "Constrained ×\nSomewhat dangerous",
                    "Constrained ×\nSomewhat safe",
                    "Constrained ×\nVery safe"))
acie_constr_safety$color <- ifelse(acie_constr_safety$estimate == 0, 0, 1)

p1 <- ggplot(acie_constr_safety, aes(x = estimate, y = term, color = factor(color))) +
  geom_point(size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0) +
  xlab("Effect on Pr(preferring to be born\nand grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Institutional Checks vs. Public Safety") +
  coord_cartesian(xlim = c(-.25, .25))
p1 <- p1 + scale_color_manual(values = c("grey60", "black"))
p1 <- p1 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 11, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white")
  )

### Analyze interaction effects between institutional checks and respondent wealth (Figure 5B) ----
## Create a new interaction variable
df_cj$constr_wealth <- interaction(df_cj$`Institutional Checks`, df_cj$`Respondent Wealth`, sep = " × ")
df_cj$constr_wealth <- relevel(df_cj$constr_wealth, ref = "Constrained × Poor")

## Estimate the average component interaction effects
acie_constr_wealth_est <- 
  lm_robust(selected ~ constr_wealth + `Leader Selection` + `Civil Liberties` +
              `Corruption in Politics` + `National Economy` + `Public Safety` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, clusters = id, fixed_effects = country)
acie_constr_wealth <- acie_constr_wealth_est %>% 
  tidy() %>% 
  add_row(term = "Constrained × Poor", estimate = 0) %>% 
  slice(1:8, 27)
acie_constr_wealth$term <- 
  factor(acie_constr_wealth$term,
         levels = c("constr_wealthUnconstrained × Poor",
                    "constr_wealthUnconstrained × Average",
                    "constr_wealthUnconstrained × Wealthy",
                    "constr_wealthPartially constrained × Poor",
                    "constr_wealthPartially constrained × Average",
                    "constr_wealthPartially constrained × Wealthy",
                    "Constrained × Poor",
                    "constr_wealthConstrained × Average",
                    "constr_wealthConstrained × Wealthy"),
         labels = c("Unconstrained ×\nPoor",
                    "Unconstrained ×\nAverage",
                    "Unconstrained ×\nWealthy",
                    "Partially constrained ×\nPoor",
                    "Partially constrained ×\nAverage",
                    "Partially constrained ×\nWealthy",
                    "Constrained ×\nPoor",
                    "Constrained ×\nAverage",
                    "Constrained ×\nWealthy"))
acie_constr_wealth$color <- ifelse(acie_constr_wealth$estimate == 0, 0, 1)

p2 <- ggplot(acie_constr_wealth, aes(x = estimate, y = term, color = factor(color))) +
  geom_point(size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0) +
  xlab("Effect on Pr(preferring to be born\nand grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Institutional Checks vs. Respondent Wealth") +
  coord_cartesian(xlim = c(-.25, .25))
p2 <- p2 + scale_color_manual(values = c("grey60", "black"))
p2 <- p2 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 11, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white")
  )

plot_grid(p1, p2, 
          labels = "AUTO", ncol = 2, 
          label_fontfamily = "Times")
ggsave("Figure 5.pdf", width = 10, height = 6)

## Present results in table form (Tables S16-S17)
acies_safety <- subset(acie_constr_safety, select = c("term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
acies_safety$N <- 36900
kable(acies_safety, digits = 3, "latex")

acies_wealth <- subset(acie_constr_wealth, select = c("term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
acies_wealth$N <- 36900
kable(acies_wealth, digits = 3, "latex")

### Analyze interaction effects between institutional checks and national economy (Figure S13) ----
## Create a new interaction variable
df_cj$constr_econ <- interaction(df_cj$`Institutional Checks`, df_cj$`National Economy`, sep = " × ")
df_cj$constr_econ <- relevel(df_cj$constr_econ, ref = "Constrained × Low income")

## Estimate the average component interaction effects
acie_constr_econ_est <- 
  lm_robust(selected ~ constr_econ + `Leader Selection` + `Civil Liberties` +
              `Corruption in Politics` + `Respondent Wealth` + `Public Safety` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, clusters = id, fixed_effects = country)
acie_constr_econ <- acie_constr_econ_est %>% 
  tidy() %>% 
  add_row(term = "Constrained × Low income", estimate = 0) %>% 
  slice(1:8, 27)
acie_constr_econ$term <- 
  factor(acie_constr_econ$term,
         levels = c("constr_econUnconstrained × Low income",
                    "constr_econUnconstrained × Middle income",
                    "constr_econUnconstrained × High income",
                    "constr_econPartially constrained × Low income",
                    "constr_econPartially constrained × Middle income",
                    "constr_econPartially constrained × High income",
                    "Constrained × Low income",
                    "constr_econConstrained × Middle income",
                    "constr_econConstrained × High income"),
         labels = c("Unconstrained ×\nLow income",
                    "Unconstrained ×\nMiddle income",
                    "Unconstrained ×\nHigh income",
                    "Partially constrained ×\nLow income",
                    "Partially constrained ×\nMiddle income",
                    "Partially constrained ×\nHigh income",
                    "Constrained ×\nLow income",
                    "Constrained ×\nMiddle income",
                    "Constrained ×\nHigh income"))
acie_constr_econ$color <- ifelse(acie_constr_econ$estimate == 0, 0, 1)

p3 <- ggplot(acie_constr_econ, aes(x = estimate, y = term, color = factor(color))) +
  geom_point(size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0) +
  xlab("Effect on Pr(preferring to be born\nand grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Institutional Checks vs. National Economy") +
  coord_cartesian(xlim = c(-.25, .25))
p3 <- p3 + scale_color_manual(values = c("grey60", "black"))
p3 <- p3 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 11, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white")
  )
p3
ggsave("Figure S13.pdf", width = 6, height = 6)

## Present results in table form (Table S34)
acies_econ <- subset(acie_constr_econ, select = c("term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
acies_econ$N <- 36900
kable(acies_econ, digits = 3, "latex")

### Analyze interaction effects between institutional checks and health care (Figure S16) ----
## Create a new interaction variable
df_cj$constr_health <- interaction(df_cj$`Institutional Checks`, df_cj$`Health Care`, sep = " × ")
df_cj$constr_health <- relevel(df_cj$constr_health, ref = "Constrained × For the privileged")

## Estimate the average component interaction effects
acie_constr_health_est <- 
  lm_robust(selected ~ constr_health + `Leader Selection` + `Civil Liberties` +
              `Corruption in Politics` + `Respondent Wealth` + `National Economy` + 
              `Public Safety` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, clusters = id, fixed_effects = country)
acie_constr_health <- acie_constr_health_est %>% 
  tidy() %>% 
  add_row(term = "Constrained × For the privileged", estimate = 0) %>% 
  slice(1:5, 25)
acie_constr_health$term <- 
  factor(acie_constr_health$term,
         levels = c("constr_healthUnconstrained × For the privileged",
                    "constr_healthUnconstrained × Mostly accessible",
                    "constr_healthPartially constrained × For the privileged",
                    "constr_healthPartially constrained × Mostly accessible",
                    "Constrained × For the privileged",
                    "constr_healthConstrained × Mostly accessible"),
         labels = c("Unconstrained ×\nHealth care for the privileged",
                    "Unconstrained ×\nHealth care mostly accessible",
                    "Partially constrained ×\nHealth care for the privileged",
                    "Partially constrained ×\nHealth care mostly accessible",
                    "Constrained ×\nHealth care for the privileged",
                    "Constrained ×\nHealth care mostly accessible"))
acie_constr_health$color <- ifelse(acie_constr_health$estimate == 0, 0, 1)

p4 <- ggplot(acie_constr_health, aes(x = estimate, y = term, color = factor(color))) +
  geom_point(size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0) +
  xlab("Effect on Pr(preferring to be born\nand grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Institutional Checks vs. Health Care") +
  coord_cartesian(xlim = c(-.25, .25))
p4 <- p4 + scale_color_manual(values = c("grey60", "black"))
p4 <- p4 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 11, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white")
  )
p4
ggsave("Figure S16.pdf", width = 6, height = 6)

## Present results in table form (Table S37)
acies_health <- subset(acie_constr_health, select = c("term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
acies_health$N <- 36900
kable(acies_health, digits = 3, "latex")

### Analyze the distribution of importance of democracy (Figure S24) ----
df_temp <- df_cj %>% slice(which(row_number() %% 6 == 1))
summary(df_temp$democracy_impt_1)
sd(df_temp$democracy_impt_1)

p <- ggplot(data = df_temp, aes(x = democracy_impt_1)) +
  geom_histogram(binwidth = 1, color = "black", fill = "gray80") + 
  scale_x_continuous(n.breaks = 9) + 
  coord_cartesian(xlim = c(0.75, 10.25)) +
  scale_color_grey() +
  xlab("Importance of living in a democracy\n(1 = not at all important; 10 = absolutely important)") +
  ylab("Number of respondents")
p <- p + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 10, colour = "black"),
    axis.text.y = element_text(size = 10, colour = "black"),
    axis.ticks = element_line(colour = "grey50"),
    strip.text.x = element_text(size = 12)
  )
p
ggsave("Figure S24.pdf", width = 8, height = 4)

### Analyze interaction effects between leader selection and public safety by importance of democracy (Figure S25) ----
## Importance of democracy above mean (8.86) vs. below mean
df_cj$demo_impt_bin_2 <- ifelse(df_cj$democracy_impt_1 > mean(df_cj$democracy_impt_1), "Importance of Democracy Above Mean", "Below")
df_cj$demo_impt_bin_2 <- factor(df_cj$demo_impt_bin_2, levels = c("Below", "Importance of Democracy Above Mean"))

## Estimate the average component interaction effects among respondents whose importance of democracy is above mean
acie_election_safety_est_1 <- 
  lm_robust(selected ~ election_safety + `Civil Liberties` + `Institutional Checks` +
              `Corruption in Politics` + `National Economy` + `Respondent Wealth` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, subset = demo_impt_bin_2 == "Importance of Democracy Above Mean",
            clusters = id, fixed_effects = country)
acie_election_safety_1 <- acie_election_safety_est_1 %>% 
  tidy() %>% 
  add_row(term = "Free and fair elections × Very dangerous", estimate = 0) %>% 
  slice(1:7, 23) %>% 
  mutate("Importance of Democracy" = "Importance of democracy above mean (9-10)")
acie_election_safety_1$term <- 
  factor(acie_election_safety_1$term,
         levels = c("election_safetyNo free and fair elections × Very dangerous",
                    "election_safetyNo free and fair elections × Somewhat dangerous",
                    "election_safetyNo free and fair elections × Somewhat safe",
                    "election_safetyNo free and fair elections × Very safe",
                    "Free and fair elections × Very dangerous",
                    "election_safetyFree and fair elections × Somewhat dangerous",
                    "election_safetyFree and fair elections × Somewhat safe",
                    "election_safetyFree and fair elections × Very safe"),
         labels = c("No free and fair elections ×\nVery dangerous",
                    "No free and fair elections ×\nSomewhat dangerous",
                    "No free and fair elections ×\nSomewhat safe",
                    "No free and fair elections ×\nVery safe",
                    "Free and fair elections ×\nVery dangerous",
                    "Free and fair elections ×\nSomewhat dangerous",
                    "Free and fair elections ×\nSomewhat safe",
                    "Free and fair elections ×\nVery safe"))
acie_election_safety_1$color <- ifelse(acie_election_safety_1$estimate == 0, 0, 1)

## Estimate the average component interaction effects among respondents whose importance of democracy is below mean
acie_election_safety_est_2 <- 
  lm_robust(selected ~ election_safety + `Civil Liberties` + `Institutional Checks` +
              `Corruption in Politics` + `National Economy` + `Respondent Wealth` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, subset = demo_impt_bin_2 == "Below",
            clusters = id, fixed_effects = country)
acie_election_safety_2 <- acie_election_safety_est_2 %>% 
  tidy() %>% 
  add_row(term = "Free and fair elections × Very dangerous", estimate = 0) %>% 
  slice(1:7, 23) %>% 
  mutate("Importance of Democracy" = "Importance of democracy below mean (1-8)")
acie_election_safety_2$term <- 
  factor(acie_election_safety_2$term,
         levels = c("election_safetyNo free and fair elections × Very dangerous",
                    "election_safetyNo free and fair elections × Somewhat dangerous",
                    "election_safetyNo free and fair elections × Somewhat safe",
                    "election_safetyNo free and fair elections × Very safe",
                    "Free and fair elections × Very dangerous",
                    "election_safetyFree and fair elections × Somewhat dangerous",
                    "election_safetyFree and fair elections × Somewhat safe",
                    "election_safetyFree and fair elections × Very safe"),
         labels = c("No free and fair elections ×\nVery dangerous",
                    "No free and fair elections ×\nSomewhat dangerous",
                    "No free and fair elections ×\nSomewhat safe",
                    "No free and fair elections ×\nVery safe",
                    "Free and fair elections ×\nVery dangerous",
                    "Free and fair elections ×\nSomewhat dangerous",
                    "Free and fair elections ×\nSomewhat safe",
                    "Free and fair elections ×\nVery safe"))
acie_election_safety_2$color <- ifelse(acie_election_safety_2$estimate == 0, 0, 1)

## Merge the estimates and plot the results
acie_election_safety_merged <- bind_rows(acie_election_safety_1, acie_election_safety_2)
p1 <- ggplot(acie_election_safety_merged, aes(x = estimate, y = term, color = factor(`Importance of Democracy`),
                                              shape = factor(`Importance of Democracy`),
                                              linetype = factor(`Importance of Democracy`))) +
  geom_point(size = 2, position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0,
                position = position_dodge(width = 0.4)) +
  scale_shape_manual(values = c(17, 16)) +
  scale_linetype_manual(values = c("solid", "solid")) +
  xlab("Effect on Pr(preferring to be born and grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Free and Fair Elections vs. Public Safety") +
  coord_cartesian(xlim = c(-.26, .26))
p1 <- p1 + scale_color_manual(values = c("grey60", "black"))
p1 <- p1 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 13, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white"),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.box.background = element_rect(color = "black"),
    legend.title = element_blank(),
    legend.key.size = unit(1.5, "line"),
    legend.key.height = unit(0, "cm")
  ) +
  guides(color = guide_legend(nrow = 1, byrow = TRUE))
p1
ggsave("Figure S25.pdf", width = 9, height = 6)

## Present results in table form (Table S53)
acie_election_safety_merged$Subgroup <- 
  ifelse(acie_election_safety_merged$`Importance of Democracy` == "Importance of democracy above mean (9-10)", 
         "Above mean", "Below mean")
acies_safety <- subset(acie_election_safety_merged, select = c("Subgroup", "term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
table(df_cj$demo_impt_bin_2)
acies_safety$N <- ifelse(acies_safety$Subgroup == "Above mean", 26364, 10536)
kable(acies_safety, digits = 3, "latex")

### Analyze interaction effects between civil liberties and public safety by importance of democracy (Figure S26) ----
## Estimate the average component interaction effects among respondents whose importance of democracy is above mean
acie_civil_safety_est_1 <- 
  lm_robust(selected ~ civil_safety + `Leader Selection` + `Institutional Checks` +
              `Corruption in Politics` + `National Economy` + `Respondent Wealth` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, subset = demo_impt_bin_2 == "Importance of Democracy Above Mean",
            clusters = id, fixed_effects = country)
acie_civil_safety_1 <- acie_civil_safety_est_1 %>% 
  tidy() %>% 
  add_row(term = "Free × Very dangerous", estimate = 0) %>% 
  slice(1:11, 29) %>% 
  mutate("Importance of Democracy" = "Importance of democracy above mean (9-10)")
acie_civil_safety_1$term <- 
  factor(acie_civil_safety_1$term,
         levels = c("civil_safetyRepressed × Very dangerous",
                    "civil_safetyRepressed × Somewhat dangerous",
                    "civil_safetyRepressed × Somewhat safe",
                    "civil_safetyRepressed × Very safe",
                    "civil_safetyPartially free × Very dangerous",
                    "civil_safetyPartially free × Somewhat dangerous",
                    "civil_safetyPartially free × Somewhat safe",
                    "civil_safetyPartially free × Very safe",
                    "Free × Very dangerous",
                    "civil_safetyFree × Somewhat dangerous",
                    "civil_safetyFree × Somewhat safe",
                    "civil_safetyFree × Very safe"),
         labels = c("Repressed ×\nVery dangerous",
                    "Repressed ×\nSomewhat dangerous",
                    "Repressed ×\nSomewhat safe",
                    "Repressed ×\nVery safe",
                    "Partially free ×\nVery dangerous",
                    "Partially free ×\nSomewhat dangerous",
                    "Partially free ×\nSomewhat safe",
                    "Partially free ×\nVery safe",
                    "Free ×\nVery dangerous",
                    "Free ×\nSomewhat dangerous",
                    "Free ×\nSomewhat safe",
                    "Free ×\nVery safe"))
acie_civil_safety_1$color <- ifelse(acie_civil_safety_1$estimate == 0, 0, 1)

## Estimate the average component interaction effects among respondents whose importance of democracy is below mean
acie_civil_safety_est_2 <- 
  lm_robust(selected ~ civil_safety + `Leader Selection` + `Institutional Checks` +
              `Corruption in Politics` + `National Economy` + `Respondent Wealth` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, subset = demo_impt_bin_2 == "Below",
            clusters = id, fixed_effects = country)
acie_civil_safety_2 <- acie_civil_safety_est_2 %>% 
  tidy() %>% 
  add_row(term = "Free × Very dangerous", estimate = 0) %>% 
  slice(1:11, 29) %>% 
  mutate("Importance of Democracy" = "Importance of democracy below mean (1-8)")
acie_civil_safety_2$term <- 
  factor(acie_civil_safety_2$term,
         levels = c("civil_safetyRepressed × Very dangerous",
                    "civil_safetyRepressed × Somewhat dangerous",
                    "civil_safetyRepressed × Somewhat safe",
                    "civil_safetyRepressed × Very safe",
                    "civil_safetyPartially free × Very dangerous",
                    "civil_safetyPartially free × Somewhat dangerous",
                    "civil_safetyPartially free × Somewhat safe",
                    "civil_safetyPartially free × Very safe",
                    "Free × Very dangerous",
                    "civil_safetyFree × Somewhat dangerous",
                    "civil_safetyFree × Somewhat safe",
                    "civil_safetyFree × Very safe"),
         labels = c("Repressed ×\nVery dangerous",
                    "Repressed ×\nSomewhat dangerous",
                    "Repressed ×\nSomewhat safe",
                    "Repressed ×\nVery safe",
                    "Partially free ×\nVery dangerous",
                    "Partially free ×\nSomewhat dangerous",
                    "Partially free ×\nSomewhat safe",
                    "Partially free ×\nVery safe",
                    "Free ×\nVery dangerous",
                    "Free ×\nSomewhat dangerous",
                    "Free ×\nSomewhat safe",
                    "Free ×\nVery safe"))
acie_civil_safety_2$color <- ifelse(acie_civil_safety_2$estimate == 0, 0, 1)

## Merge the estimates and plot the results
acie_civil_safety_merged <- bind_rows(acie_civil_safety_1, acie_civil_safety_2)
p2 <- ggplot(acie_civil_safety_merged, aes(x = estimate, y = term, color = factor(`Importance of Democracy`),
                                           shape = factor(`Importance of Democracy`),
                                           linetype = factor(`Importance of Democracy`))) +
  geom_point(size = 2, position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0,
                position = position_dodge(width = 0.4)) +
  scale_shape_manual(values = c(17, 16)) +
  scale_linetype_manual(values = c("solid", "solid")) +
  xlab("Effect on Pr(preferring to be born and grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Civil Liberties vs. Public Safety") +
  coord_cartesian(xlim = c(-.26, .26))
p2 <- p2 + scale_color_manual(values = c("grey60", "black"))
p2 <- p2 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 13, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white"),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.box.background = element_rect(color = "black"),
    legend.title = element_blank(),
    legend.key.size = unit(1.5, "line"),
    legend.key.height = unit(0, "cm")
  ) +
  guides(color = guide_legend(nrow = 1, byrow = TRUE))
p2
ggsave("Figure S26.pdf", width = 9, height = 6)

## Present results in table form (Table S54)
acie_civil_safety_merged$Subgroup <- 
  ifelse(acie_civil_safety_merged$`Importance of Democracy` == "Importance of democracy above mean (9-10)", 
         "Above mean", "Below mean")
acies_safety <- subset(acie_civil_safety_merged, select = c("Subgroup", "term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
table(df_cj$demo_impt_bin_2)
acies_safety$N <- ifelse(acies_safety$Subgroup == "Above mean", 26364, 10536)
kable(acies_safety, digits = 3, "latex")

### Analyze interaction effects between institutional checks and public safety by importance of democracy (Figure S27) ----
## Estimate the average component interaction effects among respondents whose importance of democracy is above mean
acie_constr_safety_est_1 <- 
  lm_robust(selected ~ constr_safety + `Leader Selection` + `Civil Liberties` +
              `Corruption in Politics` + `National Economy` + `Respondent Wealth` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, subset = demo_impt_bin_2 == "Importance of Democracy Above Mean",
            clusters = id, fixed_effects = country)
acie_constr_safety_1 <- acie_constr_safety_est_1 %>% 
  tidy() %>% 
  add_row(term = "Constrained × Very dangerous", estimate = 0) %>% 
  slice(1:11, 29) %>% 
  mutate("Importance of Democracy" = "Importance of democracy above mean (9-10)")
acie_constr_safety_1$term <- 
  factor(acie_constr_safety_1$term,
         levels = c("constr_safetyUnconstrained × Very dangerous",
                    "constr_safetyUnconstrained × Somewhat dangerous",
                    "constr_safetyUnconstrained × Somewhat safe",
                    "constr_safetyUnconstrained × Very safe",
                    "constr_safetyPartially constrained × Very dangerous",
                    "constr_safetyPartially constrained × Somewhat dangerous",
                    "constr_safetyPartially constrained × Somewhat safe",
                    "constr_safetyPartially constrained × Very safe",
                    "Constrained × Very dangerous",
                    "constr_safetyConstrained × Somewhat dangerous",
                    "constr_safetyConstrained × Somewhat safe",
                    "constr_safetyConstrained × Very safe"),
         labels = c("Unconstrained ×\nVery dangerous",
                    "Unconstrained ×\nSomewhat dangerous",
                    "Unconstrained ×\nSomewhat safe",
                    "Unconstrained ×\nVery safe",
                    "Partially constrained ×\nVery dangerous",
                    "Partially constrained ×\nSomewhat dangerous",
                    "Partially constrained ×\nSomewhat safe",
                    "Partially constrained ×\nVery safe",
                    "Constrained ×\nVery dangerous",
                    "Constrained ×\nSomewhat dangerous",
                    "Constrained ×\nSomewhat safe",
                    "Constrained ×\nVery safe"))
acie_constr_safety_1$color <- ifelse(acie_constr_safety_1$estimate == 0, 0, 1)

## Estimate the average component interaction effects among respondents whose importance of democracy is below mean
acie_constr_safety_est_2 <- 
  lm_robust(selected ~ constr_safety + `Leader Selection` + `Civil Liberties` +
              `Corruption in Politics` + `National Economy` + `Respondent Wealth` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, subset = demo_impt_bin_2 == "Below",
            clusters = id, fixed_effects = country)
acie_constr_safety_2 <- acie_constr_safety_est_2 %>% 
  tidy() %>% 
  add_row(term = "Constrained × Very dangerous", estimate = 0) %>% 
  slice(1:11, 29) %>% 
  mutate("Importance of Democracy" = "Importance of democracy below mean (1-8)")
acie_constr_safety_2$term <- 
  factor(acie_constr_safety_2$term,
         levels = c("constr_safetyUnconstrained × Very dangerous",
                    "constr_safetyUnconstrained × Somewhat dangerous",
                    "constr_safetyUnconstrained × Somewhat safe",
                    "constr_safetyUnconstrained × Very safe",
                    "constr_safetyPartially constrained × Very dangerous",
                    "constr_safetyPartially constrained × Somewhat dangerous",
                    "constr_safetyPartially constrained × Somewhat safe",
                    "constr_safetyPartially constrained × Very safe",
                    "Constrained × Very dangerous",
                    "constr_safetyConstrained × Somewhat dangerous",
                    "constr_safetyConstrained × Somewhat safe",
                    "constr_safetyConstrained × Very safe"),
         labels = c("Unconstrained ×\nVery dangerous",
                    "Unconstrained ×\nSomewhat dangerous",
                    "Unconstrained ×\nSomewhat safe",
                    "Unconstrained ×\nVery safe",
                    "Partially constrained ×\nVery dangerous",
                    "Partially constrained ×\nSomewhat dangerous",
                    "Partially constrained ×\nSomewhat safe",
                    "Partially constrained ×\nVery safe",
                    "Constrained ×\nVery dangerous",
                    "Constrained ×\nSomewhat dangerous",
                    "Constrained ×\nSomewhat safe",
                    "Constrained ×\nVery safe"))
acie_constr_safety_2$color <- ifelse(acie_constr_safety_2$estimate == 0, 0, 1)

## Merge the estimates and plot the results
acie_constr_safety_merged <- bind_rows(acie_constr_safety_1, acie_constr_safety_2)
p3 <- ggplot(acie_constr_safety_merged, aes(x = estimate, y = term, color = factor(`Importance of Democracy`),
                                            shape = factor(`Importance of Democracy`),
                                            linetype = factor(`Importance of Democracy`))) +
  geom_point(size = 2, position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0,
                position = position_dodge(width = 0.4)) +
  scale_shape_manual(values = c(17, 16)) +
  scale_linetype_manual(values = c("solid", "solid")) +
  xlab("Effect on Pr(preferring to be born and grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Institutional Checks vs. Public Safety") +
  coord_cartesian(xlim = c(-.26, .26))
p3 <- p3 + scale_color_manual(values = c("grey60", "black"))
p3 <- p3 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 13, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white"),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.box.background = element_rect(color = "black"),
    legend.title = element_blank(),
    legend.key.size = unit(1.5, "line"),
    legend.key.height = unit(0, "cm")
  ) +
  guides(color = guide_legend(nrow = 1, byrow = TRUE))
p3
ggsave("Figure S27.pdf", width = 9, height = 6)

## Present results in table form (Table S55)
acie_constr_safety_merged$Subgroup <- 
  ifelse(acie_constr_safety_merged$`Importance of Democracy` == "Importance of democracy above mean (9-10)", 
         "Above mean", "Below mean")
acies_safety <- subset(acie_constr_safety_merged, select = c("Subgroup", "term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
table(df_cj$demo_impt_bin_2)
acies_safety$N <- ifelse(acies_safety$Subgroup == "Above mean", 26364, 10536)
kable(acies_safety, digits = 3, "latex")

### Analyze interaction effects between leader selection and respondent wealth by importance of democracy (Figure S28) ----
## Estimate the average component interaction effects among respondents whose importance of democracy is above mean
acie_election_wealth_est_1 <- 
  lm_robust(selected ~ election_wealth + `Civil Liberties` + `Institutional Checks` +
              `Corruption in Politics` + `National Economy` + `Public Safety` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, subset = demo_impt_bin_2 == "Importance of Democracy Above Mean",
            clusters = id, fixed_effects = country)
acie_election_wealth_1 <- acie_election_wealth_est_1 %>% 
  tidy() %>% 
  add_row(term = "Free and fair elections × Poor", estimate = 0) %>% 
  slice(1:5, 22) %>% 
  mutate("Importance of Democracy" = "Importance of democracy above mean (9-10)")
acie_election_wealth_1$term <- 
  factor(acie_election_wealth_1$term,
         levels = c("election_wealthNo free and fair elections × Poor",
                    "election_wealthNo free and fair elections × Average",
                    "election_wealthNo free and fair elections × Wealthy",
                    "Free and fair elections × Poor",
                    "election_wealthFree and fair elections × Average",
                    "election_wealthFree and fair elections × Wealthy"),
         labels = c("No free and fair elections ×\nPoor",
                    "No free and fair elections ×\nAverage",
                    "No free and fair elections ×\nWealthy",
                    "Free and fair elections ×\nPoor",
                    "Free and fair elections ×\nAverage",
                    "Free and fair elections ×\nWealthy"))
acie_election_wealth_1$color <- ifelse(acie_election_wealth_1$estimate == 0, 0, 1)

## Estimate the average component interaction effects among respondents whose importance of democracy is below mean
acie_election_wealth_est_2 <- 
  lm_robust(selected ~ election_wealth + `Civil Liberties` + `Institutional Checks` +
              `Corruption in Politics` + `National Economy` + `Public Safety` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, subset = demo_impt_bin_2 == "Below",
            clusters = id, fixed_effects = country)
acie_election_wealth_2 <- acie_election_wealth_est_2 %>% 
  tidy() %>% 
  add_row(term = "Free and fair elections × Poor", estimate = 0) %>% 
  slice(1:5, 22) %>% 
  mutate("Importance of Democracy" = "Importance of democracy below mean (1-8)")
acie_election_wealth_2$term <- 
  factor(acie_election_wealth_2$term,
         levels = c("election_wealthNo free and fair elections × Poor",
                    "election_wealthNo free and fair elections × Average",
                    "election_wealthNo free and fair elections × Wealthy",
                    "Free and fair elections × Poor",
                    "election_wealthFree and fair elections × Average",
                    "election_wealthFree and fair elections × Wealthy"),
         labels = c("No free and fair elections ×\nPoor",
                    "No free and fair elections ×\nAverage",
                    "No free and fair elections ×\nWealthy",
                    "Free and fair elections ×\nPoor",
                    "Free and fair elections ×\nAverage",
                    "Free and fair elections ×\nWealthy"))
acie_election_wealth_2$color <- ifelse(acie_election_wealth_2$estimate == 0, 0, 1)

## Merge the estimates and plot the results
acie_election_wealth_merged <- bind_rows(acie_election_wealth_1, acie_election_wealth_2)
p1 <- ggplot(acie_election_wealth_merged, aes(x = estimate, y = term, color = factor(`Importance of Democracy`),
                                              shape = factor(`Importance of Democracy`),
                                              linetype = factor(`Importance of Democracy`))) +
  geom_point(size = 2, position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0,
                position = position_dodge(width = 0.4)) +
  scale_shape_manual(values = c(17, 16)) +
  scale_linetype_manual(values = c("solid", "solid")) +
  xlab("Effect on Pr(preferring to be born and grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Free and Fair Elections vs. Respondent Wealth") +
  coord_cartesian(xlim = c(-.26, .26))
p1 <- p1 + scale_color_manual(values = c("grey60", "black"))
p1 <- p1 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 13, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white"),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.box.background = element_rect(color = "black"),
    legend.title = element_blank(),
    legend.key.size = unit(1.5, "line"),
    legend.key.height = unit(0, "cm")
  ) +
  guides(color = guide_legend(nrow = 1, byrow = TRUE))
p1
ggsave("Figure S28.pdf", width = 9, height = 6)

## Present results in table form (Table S56)
acie_election_wealth_merged$Subgroup <- 
  ifelse(acie_election_wealth_merged$`Importance of Democracy` == "Importance of democracy above mean (9-10)", 
         "Above mean", "Below mean")
acies_wealth <- subset(acie_election_wealth_merged, select = c("Subgroup", "term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
table(df_cj$demo_impt_bin_2)
acies_wealth$N <- ifelse(acies_wealth$Subgroup == "Above mean", 26364, 10536)
kable(acies_wealth, digits = 3, "latex")

### Analyze interaction effects between civil liberties and respondent wealth by importance of democracy (Figure S29) ----
## Estimate the average component interaction effects among respondents whose importance of democracy is above mean
acie_civil_wealth_est_1 <- 
  lm_robust(selected ~ civil_wealth + `Leader Selection` + `Institutional Checks` +
              `Corruption in Politics` + `National Economy` + `Public Safety` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, subset = demo_impt_bin_2 == "Importance of Democracy Above Mean",
            clusters = id, fixed_effects = country)
acie_civil_wealth_1 <- acie_civil_wealth_est_1 %>% 
  tidy() %>% 
  add_row(term = "Free × Poor", estimate = 0) %>% 
  slice(1:8, 27) %>% 
  mutate("Importance of Democracy" = "Importance of democracy above mean (9-10)")
acie_civil_wealth_1$term <- 
  factor(acie_civil_wealth_1$term,
         levels = c("civil_wealthRepressed × Poor",
                    "civil_wealthRepressed × Average",
                    "civil_wealthRepressed × Wealthy",
                    "civil_wealthPartially free × Poor",
                    "civil_wealthPartially free × Average",
                    "civil_wealthPartially free × Wealthy",
                    "Free × Poor",
                    "civil_wealthFree × Average",
                    "civil_wealthFree × Wealthy"),
         labels = c("Repressed ×\nPoor",
                    "Repressed ×\nAverage",
                    "Repressed ×\nWealthy",
                    "Partially free ×\nPoor",
                    "Partially free ×\nAverage",
                    "Partially free ×\nWealthy",
                    "Free ×\nPoor",
                    "Free ×\nAverage",
                    "Free ×\nWealthy"))
acie_civil_wealth_1$color <- ifelse(acie_civil_wealth_1$estimate == 0, 0, 1)

## Estimate the average component interaction effects among respondents whose importance of democracy is below mean
acie_civil_wealth_est_2 <- 
  lm_robust(selected ~ civil_wealth + `Leader Selection` + `Institutional Checks` +
              `Corruption in Politics` + `National Economy` + `Public Safety` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, subset = demo_impt_bin_2 == "Below",
            clusters = id, fixed_effects = country)
acie_civil_wealth_2 <- acie_civil_wealth_est_2 %>% 
  tidy() %>% 
  add_row(term = "Free × Poor", estimate = 0) %>% 
  slice(1:8, 27) %>% 
  mutate("Importance of Democracy" = "Importance of democracy below mean (1-8)")
acie_civil_wealth_2$term <- 
  factor(acie_civil_wealth_2$term,
         levels = c("civil_wealthRepressed × Poor",
                    "civil_wealthRepressed × Average",
                    "civil_wealthRepressed × Wealthy",
                    "civil_wealthPartially free × Poor",
                    "civil_wealthPartially free × Average",
                    "civil_wealthPartially free × Wealthy",
                    "Free × Poor",
                    "civil_wealthFree × Average",
                    "civil_wealthFree × Wealthy"),
         labels = c("Repressed ×\nPoor",
                    "Repressed ×\nAverage",
                    "Repressed ×\nWealthy",
                    "Partially free ×\nPoor",
                    "Partially free ×\nAverage",
                    "Partially free ×\nWealthy",
                    "Free ×\nPoor",
                    "Free ×\nAverage",
                    "Free ×\nWealthy"))
acie_civil_wealth_2$color <- ifelse(acie_civil_wealth_2$estimate == 0, 0, 1)

## Merge the estimates and plot the results
acie_civil_wealth_merged <- bind_rows(acie_civil_wealth_1, acie_civil_wealth_2)
p2 <- ggplot(acie_civil_wealth_merged, aes(x = estimate, y = term, color = factor(`Importance of Democracy`),
                                           shape = factor(`Importance of Democracy`),
                                           linetype = factor(`Importance of Democracy`))) +
  geom_point(size = 2, position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0,
                position = position_dodge(width = 0.4)) +
  scale_shape_manual(values = c(17, 16)) +
  scale_linetype_manual(values = c("solid", "solid")) +
  xlab("Effect on Pr(preferring to be born and grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Civil Liberties vs. Respondent Wealth") +
  coord_cartesian(xlim = c(-.26, .26))
p2 <- p2 + scale_color_manual(values = c("grey60", "black"))
p2 <- p2 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 13, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white"),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.box.background = element_rect(color = "black"),
    legend.title = element_blank(),
    legend.key.size = unit(1.5, "line"),
    legend.key.height = unit(0, "cm")
  ) +
  guides(color = guide_legend(nrow = 1, byrow = TRUE))
p2
ggsave("Figure S29.pdf", width = 9, height = 6)

## Present results in table form (Table S57)
acie_civil_wealth_merged$Subgroup <- 
  ifelse(acie_civil_wealth_merged$`Importance of Democracy` == "Importance of democracy above mean (9-10)", 
         "Above mean", "Below mean")
acies_wealth <- subset(acie_civil_wealth_merged, select = c("Subgroup", "term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
table(df_cj$demo_impt_bin_2)
acies_wealth$N <- ifelse(acies_wealth$Subgroup == "Above mean", 26364, 10536)
kable(acies_wealth, digits = 3, "latex")

### Analyze interaction effects between institutional checks and respondent wealth by importance of democracy (Figure S30) ----
## Estimate the average component interaction effects among respondents whose importance of democracy is above mean
acie_constr_wealth_est_1 <- 
  lm_robust(selected ~ constr_wealth + `Leader Selection` + `Civil Liberties` +
              `Corruption in Politics` + `National Economy` + `Public Safety` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, subset = demo_impt_bin_2 == "Importance of Democracy Above Mean",
            clusters = id, fixed_effects = country)
acie_constr_wealth_1 <- acie_constr_wealth_est_1 %>% 
  tidy() %>% 
  add_row(term = "Constrained × Poor", estimate = 0) %>% 
  slice(1:8, 27) %>% 
  mutate("Importance of Democracy" = "Importance of democracy above mean (9-10)")
acie_constr_wealth_1$term <- 
  factor(acie_constr_wealth_1$term,
         levels = c("constr_wealthUnconstrained × Poor",
                    "constr_wealthUnconstrained × Average",
                    "constr_wealthUnconstrained × Wealthy",
                    "constr_wealthPartially constrained × Poor",
                    "constr_wealthPartially constrained × Average",
                    "constr_wealthPartially constrained × Wealthy",
                    "Constrained × Poor",
                    "constr_wealthConstrained × Average",
                    "constr_wealthConstrained × Wealthy"),
         labels = c("Unconstrained ×\nPoor",
                    "Unconstrained ×\nAverage",
                    "Unconstrained ×\nWealthy",
                    "Partially constrained ×\nPoor",
                    "Partially constrained ×\nAverage",
                    "Partially constrained ×\nWealthy",
                    "Constrained ×\nPoor",
                    "Constrained ×\nAverage",
                    "Constrained ×\nWealthy"))
acie_constr_wealth_1$color <- ifelse(acie_constr_wealth_1$estimate == 0, 0, 1)

## Estimate the average component interaction effects among respondents whose importance of democracy is below mean
acie_constr_wealth_est_2 <- 
  lm_robust(selected ~ constr_wealth + `Leader Selection` + `Civil Liberties` +
              `Corruption in Politics` + `National Economy` + `Public Safety` + 
              `Health Care` + `Minority Treatment` + `Respondent Identity`, 
            data = df_cj, subset = demo_impt_bin_2 == "Below",
            clusters = id, fixed_effects = country)
acie_constr_wealth_2 <- acie_constr_wealth_est_2 %>% 
  tidy() %>% 
  add_row(term = "Constrained × Poor", estimate = 0) %>% 
  slice(1:8, 27) %>% 
  mutate("Importance of Democracy" = "Importance of democracy below mean (1-8)")
acie_constr_wealth_2$term <- 
  factor(acie_constr_wealth_2$term,
         levels = c("constr_wealthUnconstrained × Poor",
                    "constr_wealthUnconstrained × Average",
                    "constr_wealthUnconstrained × Wealthy",
                    "constr_wealthPartially constrained × Poor",
                    "constr_wealthPartially constrained × Average",
                    "constr_wealthPartially constrained × Wealthy",
                    "Constrained × Poor",
                    "constr_wealthConstrained × Average",
                    "constr_wealthConstrained × Wealthy"),
         labels = c("Unconstrained ×\nPoor",
                    "Unconstrained ×\nAverage",
                    "Unconstrained ×\nWealthy",
                    "Partially constrained ×\nPoor",
                    "Partially constrained ×\nAverage",
                    "Partially constrained ×\nWealthy",
                    "Constrained ×\nPoor",
                    "Constrained ×\nAverage",
                    "Constrained ×\nWealthy"))
acie_constr_wealth_2$color <- ifelse(acie_constr_wealth_2$estimate == 0, 0, 1)

## Merge the estimates and plot the results
acie_constr_wealth_merged <- bind_rows(acie_constr_wealth_1, acie_constr_wealth_2)
p3 <- ggplot(acie_constr_wealth_merged, aes(x = estimate, y = term, color = factor(`Importance of Democracy`),
                                            shape = factor(`Importance of Democracy`),
                                            linetype = factor(`Importance of Democracy`))) +
  geom_point(size = 2, position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, width = 0,
                position = position_dodge(width = 0.4)) +
  scale_shape_manual(values = c(17, 16)) +
  scale_linetype_manual(values = c("solid", "solid")) +
  xlab("Effect on Pr(preferring to be born and grow up in the country)") +
  ylab("") +
  geom_vline(xintercept = 0, color = "grey") +
  ggtitle("Institutional Checks vs. Respondent Wealth") +
  coord_cartesian(xlim = c(-.26, .26))
p3 <- p3 + scale_color_manual(values = c("grey60", "black"))
p3 <- p3 + theme_bw(base_size = 12, base_family = "Times") %+replace%
  theme(
    plot.title = element_text(size = 13, hjust = .5, face = "bold.italic"),
    axis.text.x = element_text(size = 11, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 11, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.y = element_text(size = 12, angle = 90, vjust = .01, hjust = .1),
    strip.text.x = element_text(size = 12),
    strip.background = element_rect(colour = "black", fill = "white"),
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.background = element_blank(),
    legend.box.background = element_rect(color = "black"),
    legend.title = element_blank(),
    legend.key.size = unit(1.5, "line"),
    legend.key.height = unit(0, "cm")
  ) +
  guides(color = guide_legend(nrow = 1, byrow = TRUE))
p3
ggsave("Figure S30.pdf", width = 9, height = 6)

## Present results in table form (Table S58)
acie_constr_wealth_merged$Subgroup <- 
  ifelse(acie_constr_wealth_merged$`Importance of Democracy` == "Importance of democracy above mean (9-10)", 
         "Above mean", "Below mean")
acies_wealth <- subset(acie_constr_wealth_merged, select = c("Subgroup", "term", "estimate", "std.error", "p.value", "conf.low", "conf.high"))
table(df_cj$demo_impt_bin_2)
acies_wealth$N <- ifelse(acies_wealth$Subgroup == "Above mean", 26364, 10536)
kable(acies_wealth, digits = 3, "latex")

### Analyze marginal means of interaction terms (Figures S17-S22) ----
## Reverse the factors
df_cj$`Leader Selection` <- fct_rev(df_cj$`Leader Selection`)
df_cj$`Civil Liberties` <- fct_rev(df_cj$`Civil Liberties`)
df_cj$`Institutional Checks` <- fct_rev(df_cj$`Institutional Checks`)
df_cj$`National Economy` <- fct_rev(df_cj$`National Economy`)
df_cj$`Respondent Wealth` <- fct_rev(df_cj$`Respondent Wealth`)
df_cj$`Corruption in Politics` <- fct_rev(df_cj$`Corruption in Politics`)
df_cj$`Health Care` <- fct_rev(df_cj$`Health Care`)
df_cj$`Public Safety` <- fct_rev(df_cj$`Public Safety`)
df_cj$`Minority Treatment` <- fct_rev(df_cj$`Minority Treatment`)
df_cj$`Respondent Identity` <- fct_rev(df_cj$`Respondent Identity`)

## Leader selection vs. public safety (Figure S17)
# Create a new interaction variable
df_cj$election_safety <- interaction(df_cj$`Leader Selection`, df_cj$`Public Safety`, sep = " × ")

# Estimate the MMs
mm_interact_EG <- 
  mm(data = subset(df_cj, country == "Egypt"), 
     formula = selected ~ election_safety,
     id = ~id) %>% 
  mutate(BY = "Egypt")

mm_interact_IN <- 
  mm(data = subset(df_cj, country == "India"), 
     formula = selected ~ election_safety,
     id = ~id) %>% 
  mutate(BY = "India")

mm_interact_IT <- 
  mm(data = subset(df_cj, country == "Italy"), 
     formula = selected ~ election_safety,
     id = ~id) %>% 
  mutate(BY = "Italy")

mm_interact_JP <- 
  mm(data = subset(df_cj, country == "Japan"), 
     formula = selected ~ election_safety,
     id = ~id) %>% 
  mutate(BY = "Japan")

mm_interact_TH <- 
  mm(data = subset(df_cj, country == "Thailand"), 
     formula = selected ~ election_safety,
     id = ~id) %>% 
  mutate(BY = "Thailand")

mm_interact_US <- 
  mm(data = subset(df_cj, country == "United States"), 
     formula = selected ~ election_safety,
     id = ~id) %>% 
  mutate(BY = "United States")

mm_interact <- bind_rows(mm_interact_EG, mm_interact_IN, mm_interact_IT, 
                         mm_interact_JP, mm_interact_TH, mm_interact_US)

# Visualize the MMs
p <- ggplot(mm_interact, aes(x = estimate, y = level)) +
  geom_point(size = 1.5) +
  geom_errorbar(aes(xmin = lower, xmax = upper), linewidth = .5, width = 0) +
  ggplot2::facet_wrap(~BY, ncol = 3) +
  xlab("Marginal mean") +
  ylab("") +
  geom_vline(xintercept = 0.5, color = "grey")
p <- p + theme_bw(base_size = 8, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 8, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 8, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 8)
  )
p <- p + coord_cartesian(xlim = c(0.2, 0.8))
p
ggsave("Figure S17.pdf", width = 6.5, height = 8.5)

## Present results in table form (Tables S38-S40)
mm_interact <- subset(mm_interact, select = c("BY", "level", "estimate", "std.error", "lower", "upper"))
mm_interact <- mm_interact %>% mutate(N = case_when(
  BY == "Egypt" ~ nrow(df_EG),
  BY == "India" ~ nrow(df_IN),
  BY == "Italy" ~ nrow(df_IT),
  BY == "Japan" ~ nrow(df_JP),
  BY == "Thailand" ~ nrow(df_TH),
  BY == "United States" ~ nrow(df_US)
))
kable(mm_interact, digits = 3, "latex", row.names = F)

## Leader selection vs. respondent wealth (Figure S18)
# Create a new interaction variable
df_cj$election_wealth <- interaction(df_cj$`Leader Selection`, df_cj$`Respondent Wealth`, sep = " × ")

# Estimate the MMs
mm_interact_EG <- 
  mm(data = subset(df_cj, country == "Egypt"), 
     formula = selected ~ election_wealth,
     id = ~id) %>% 
  mutate(BY = "Egypt")

mm_interact_IN <- 
  mm(data = subset(df_cj, country == "India"), 
     formula = selected ~ election_wealth,
     id = ~id) %>% 
  mutate(BY = "India")

mm_interact_IT <- 
  mm(data = subset(df_cj, country == "Italy"), 
     formula = selected ~ election_wealth,
     id = ~id) %>% 
  mutate(BY = "Italy")

mm_interact_JP <- 
  mm(data = subset(df_cj, country == "Japan"), 
     formula = selected ~ election_wealth,
     id = ~id) %>% 
  mutate(BY = "Japan")

mm_interact_TH <- 
  mm(data = subset(df_cj, country == "Thailand"), 
     formula = selected ~ election_wealth,
     id = ~id) %>% 
  mutate(BY = "Thailand")

mm_interact_US <- 
  mm(data = subset(df_cj, country == "United States"), 
     formula = selected ~ election_wealth,
     id = ~id) %>% 
  mutate(BY = "United States")

mm_interact <- bind_rows(mm_interact_EG, mm_interact_IN, mm_interact_IT, 
                         mm_interact_JP, mm_interact_TH, mm_interact_US)

# Visualize the MMs
p <- ggplot(mm_interact, aes(x = estimate, y = level)) +
  geom_point(size = 1.5) +
  geom_errorbar(aes(xmin = lower, xmax = upper), linewidth = .5, width = 0) +
  ggplot2::facet_wrap(~BY, ncol = 3) +
  xlab("Marginal mean") +
  ylab("") +
  geom_vline(xintercept = 0.5, color = "grey")
p <- p + theme_bw(base_size = 8, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 8, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 8, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 8)
  )
p <- p + coord_cartesian(xlim = c(0.2, 0.8))
p
ggsave("Figure S18.pdf", width = 6.5, height = 8.5)

## Present results in table form (Tables S41-S43)
mm_interact <- subset(mm_interact, select = c("BY", "level", "estimate", "std.error", "lower", "upper"))
mm_interact <- mm_interact %>% mutate(N = case_when(
  BY == "Egypt" ~ nrow(df_EG),
  BY == "India" ~ nrow(df_IN),
  BY == "Italy" ~ nrow(df_IT),
  BY == "Japan" ~ nrow(df_JP),
  BY == "Thailand" ~ nrow(df_TH),
  BY == "United States" ~ nrow(df_US)
))
kable(mm_interact, digits = 3, "latex", row.names = F)

## Civil liberties vs. public safety (Figure S19)
# Create a new interaction variable
df_cj$civil_safety <- interaction(df_cj$`Civil Liberties`, df_cj$`Public Safety`, sep = " × ")

# Estimate the MMs
mm_interact_EG <- 
  mm(data = subset(df_cj, country == "Egypt"), 
     formula = selected ~ civil_safety,
     id = ~id) %>% 
  mutate(BY = "Egypt")

mm_interact_IN <- 
  mm(data = subset(df_cj, country == "India"), 
     formula = selected ~ civil_safety,
     id = ~id) %>% 
  mutate(BY = "India")

mm_interact_IT <- 
  mm(data = subset(df_cj, country == "Italy"), 
     formula = selected ~ civil_safety,
     id = ~id) %>% 
  mutate(BY = "Italy")

mm_interact_JP <- 
  mm(data = subset(df_cj, country == "Japan"), 
     formula = selected ~ civil_safety,
     id = ~id) %>% 
  mutate(BY = "Japan")

mm_interact_TH <- 
  mm(data = subset(df_cj, country == "Thailand"), 
     formula = selected ~ civil_safety,
     id = ~id) %>% 
  mutate(BY = "Thailand")

mm_interact_US <- 
  mm(data = subset(df_cj, country == "United States"), 
     formula = selected ~ civil_safety,
     id = ~id) %>% 
  mutate(BY = "United States")

mm_interact <- bind_rows(mm_interact_EG, mm_interact_IN, mm_interact_IT, 
                         mm_interact_JP, mm_interact_TH, mm_interact_US)

# Visualize the MMs
p <- ggplot(mm_interact, aes(x = estimate, y = level)) +
  geom_point(size = 1.5) +
  geom_errorbar(aes(xmin = lower, xmax = upper), linewidth = .5, width = 0) +
  ggplot2::facet_wrap(~BY, ncol = 3) +
  xlab("Marginal mean") +
  ylab("") +
  geom_vline(xintercept = 0.5, color = "grey")
p <- p + theme_bw(base_size = 8, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 8, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 8, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 8)
  )
p <- p + coord_cartesian(xlim = c(0.2, 0.8))
p
ggsave("Figure S19.pdf", width = 6.5, height = 8.5)

## Present results in table form (Tables S44-S45)
mm_interact <- subset(mm_interact, select = c("BY", "level", "estimate", "std.error", "lower", "upper"))
mm_interact <- mm_interact %>% mutate(N = case_when(
  BY == "Egypt" ~ nrow(df_EG),
  BY == "India" ~ nrow(df_IN),
  BY == "Italy" ~ nrow(df_IT),
  BY == "Japan" ~ nrow(df_JP),
  BY == "Thailand" ~ nrow(df_TH),
  BY == "United States" ~ nrow(df_US)
))
kable(mm_interact, digits = 3, "latex", row.names = F)

## Civil liberties vs. respondent wealth (Figure S20)
# Create a new interaction variable
df_cj$civil_wealth <- interaction(df_cj$`Civil Liberties`, df_cj$`Respondent Wealth`, sep = " × ")

# Estimate the MMs
mm_interact_EG <- 
  mm(data = subset(df_cj, country == "Egypt"), 
     formula = selected ~ civil_wealth,
     id = ~id) %>% 
  mutate(BY = "Egypt")

mm_interact_IN <- 
  mm(data = subset(df_cj, country == "India"), 
     formula = selected ~ civil_wealth,
     id = ~id) %>% 
  mutate(BY = "India")

mm_interact_IT <- 
  mm(data = subset(df_cj, country == "Italy"), 
     formula = selected ~ civil_wealth,
     id = ~id) %>% 
  mutate(BY = "Italy")

mm_interact_JP <- 
  mm(data = subset(df_cj, country == "Japan"), 
     formula = selected ~ civil_wealth,
     id = ~id) %>% 
  mutate(BY = "Japan")

mm_interact_TH <- 
  mm(data = subset(df_cj, country == "Thailand"), 
     formula = selected ~ civil_wealth,
     id = ~id) %>% 
  mutate(BY = "Thailand")

mm_interact_US <- 
  mm(data = subset(df_cj, country == "United States"), 
     formula = selected ~ civil_wealth,
     id = ~id) %>% 
  mutate(BY = "United States")

mm_interact <- bind_rows(mm_interact_EG, mm_interact_IN, mm_interact_IT, 
                         mm_interact_JP, mm_interact_TH, mm_interact_US)

# Visualize the MMs
p <- ggplot(mm_interact, aes(x = estimate, y = level)) +
  geom_point(size = 1.5) +
  geom_errorbar(aes(xmin = lower, xmax = upper), linewidth = .5, width = 0) +
  ggplot2::facet_wrap(~BY, ncol = 3) +
  xlab("Marginal mean") +
  ylab("") +
  geom_vline(xintercept = 0.5, color = "grey")
p <- p + theme_bw(base_size = 8, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 8, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 8, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 8)
  )
p <- p + coord_cartesian(xlim = c(0.2, 0.8))
p
ggsave("Figure S20.pdf", width = 6.5, height = 8.5)

## Present results in table form (Tables S46-S47)
mm_interact <- subset(mm_interact, select = c("BY", "level", "estimate", "std.error", "lower", "upper"))
mm_interact <- mm_interact %>% mutate(N = case_when(
  BY == "Egypt" ~ nrow(df_EG),
  BY == "India" ~ nrow(df_IN),
  BY == "Italy" ~ nrow(df_IT),
  BY == "Japan" ~ nrow(df_JP),
  BY == "Thailand" ~ nrow(df_TH),
  BY == "United States" ~ nrow(df_US)
))
kable(mm_interact, digits = 3, "latex", row.names = F)

## Institutional checks vs. public safety (Figure S21)
# Create a new interaction variable
df_cj$constr_safety <- interaction(df_cj$`Institutional Checks`, df_cj$`Public Safety`, sep = " × ")

# Estimate the MMs
mm_interact_EG <- 
  mm(data = subset(df_cj, country == "Egypt"), 
     formula = selected ~ constr_safety,
     id = ~id) %>% 
  mutate(BY = "Egypt")

mm_interact_IN <- 
  mm(data = subset(df_cj, country == "India"), 
     formula = selected ~ constr_safety,
     id = ~id) %>% 
  mutate(BY = "India")

mm_interact_IT <- 
  mm(data = subset(df_cj, country == "Italy"), 
     formula = selected ~ constr_safety,
     id = ~id) %>% 
  mutate(BY = "Italy")

mm_interact_JP <- 
  mm(data = subset(df_cj, country == "Japan"), 
     formula = selected ~ constr_safety,
     id = ~id) %>% 
  mutate(BY = "Japan")

mm_interact_TH <- 
  mm(data = subset(df_cj, country == "Thailand"), 
     formula = selected ~ constr_safety,
     id = ~id) %>% 
  mutate(BY = "Thailand")

mm_interact_US <- 
  mm(data = subset(df_cj, country == "United States"), 
     formula = selected ~ constr_safety,
     id = ~id) %>% 
  mutate(BY = "United States")

mm_interact <- bind_rows(mm_interact_EG, mm_interact_IN, mm_interact_IT, 
                         mm_interact_JP, mm_interact_TH, mm_interact_US)

# Visualize the MMs
p <- ggplot(mm_interact, aes(x = estimate, y = level)) +
  geom_point(size = 1.5) +
  geom_errorbar(aes(xmin = lower, xmax = upper), linewidth = .5, width = 0) +
  ggplot2::facet_wrap(~BY, ncol = 3) +
  xlab("Marginal mean") +
  ylab("") +
  geom_vline(xintercept = 0.5, color = "grey")
p <- p + theme_bw(base_size = 8, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 8, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 8, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 8)
  )
p <- p + coord_cartesian(xlim = c(0.2, 0.8))
p
ggsave("Figure S21.pdf", width = 6.5, height = 8.5)

## Present results in table form (Tables S48-S49)
mm_interact <- subset(mm_interact, select = c("BY", "level", "estimate", "std.error", "lower", "upper"))
mm_interact <- mm_interact %>% mutate(N = case_when(
  BY == "Egypt" ~ nrow(df_EG),
  BY == "India" ~ nrow(df_IN),
  BY == "Italy" ~ nrow(df_IT),
  BY == "Japan" ~ nrow(df_JP),
  BY == "Thailand" ~ nrow(df_TH),
  BY == "United States" ~ nrow(df_US)
))
kable(mm_interact, digits = 3, "latex", row.names = F)

## Institutional checks vs. respondent wealth (Figure S22)
# Create a new interaction variable
df_cj$constr_wealth <- interaction(df_cj$`Institutional Checks`, df_cj$`Respondent Wealth`, sep = " × ")

# Estimate the MMs
mm_interact_EG <- 
  mm(data = subset(df_cj, country == "Egypt"), 
     formula = selected ~ constr_wealth,
     id = ~id) %>% 
  mutate(BY = "Egypt")

mm_interact_IN <- 
  mm(data = subset(df_cj, country == "India"), 
     formula = selected ~ constr_wealth,
     id = ~id) %>% 
  mutate(BY = "India")

mm_interact_IT <- 
  mm(data = subset(df_cj, country == "Italy"), 
     formula = selected ~ constr_wealth,
     id = ~id) %>% 
  mutate(BY = "Italy")

mm_interact_JP <- 
  mm(data = subset(df_cj, country == "Japan"), 
     formula = selected ~ constr_wealth,
     id = ~id) %>% 
  mutate(BY = "Japan")

mm_interact_TH <- 
  mm(data = subset(df_cj, country == "Thailand"), 
     formula = selected ~ constr_wealth,
     id = ~id) %>% 
  mutate(BY = "Thailand")

mm_interact_US <- 
  mm(data = subset(df_cj, country == "United States"), 
     formula = selected ~ constr_wealth,
     id = ~id) %>% 
  mutate(BY = "United States")

mm_interact <- bind_rows(mm_interact_EG, mm_interact_IN, mm_interact_IT, 
                         mm_interact_JP, mm_interact_TH, mm_interact_US)

# Visualize the MMs
p <- ggplot(mm_interact, aes(x = estimate, y = level)) +
  geom_point(size = 1.5) +
  geom_errorbar(aes(xmin = lower, xmax = upper), linewidth = .5, width = 0) +
  ggplot2::facet_wrap(~BY, ncol = 3) +
  xlab("Marginal mean") +
  ylab("") +
  geom_vline(xintercept = 0.5, color = "grey")
p <- p + theme_bw(base_size = 8, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 8, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 8, colour = "black", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    axis.title.x = element_text(size = 10),
    axis.title.y = element_text(size = 10, angle = 90, vjust = .01, hjust = .1),
    legend.position = "none",
    strip.text.x = element_text(size = 8)
  )
p <- p + coord_cartesian(xlim = c(0.2, 0.8))
p
ggsave("Figure S22.pdf", width = 6.5, height = 8.5)

## Present results in table form (Tables S50-S51)
mm_interact <- subset(mm_interact, select = c("BY", "level", "estimate", "std.error", "lower", "upper"))
mm_interact <- mm_interact %>% mutate(N = case_when(
  BY == "Egypt" ~ nrow(df_EG),
  BY == "India" ~ nrow(df_IN),
  BY == "Italy" ~ nrow(df_IT),
  BY == "Japan" ~ nrow(df_JP),
  BY == "Thailand" ~ nrow(df_TH),
  BY == "United States" ~ nrow(df_US)
))
kable(mm_interact, digits = 3, "latex", row.names = F)

### Attribute salience with aggregated sample (Figure S7) ----
## Reshape data
df_salience <- df_cj %>% 
  select(1:14) %>%
  mutate_if(is.factor, as.character) %>% 
  gather("attribute", "level", 4:13)

## Rename attributes
df_salience$attribute <- 
  case_match(df_salience$attribute, 
             "Leader.Selection" ~ "Leader\nSelection",
             "Civil.Liberties" ~ "Civil\nLiberties",
             "Leader.Constraints" ~ "Institutional\nChecks",
             "National.Economy" ~ "National\nEconomy",
             "Respondent.Wealth" ~ "Respondent\nWealth",
             "Public.Safety" ~ "Public\nSafety",
             "Corruption.in.Politics" ~ "Corruption\nin Politics",
             "Health.Care" ~ "Health\nCare",
             "Minority.Treatment" ~ "Minority\nTreatment",
             "Respondent.Identity" ~ "Respondent\nIdentity")

## Calculate salience scores
salience <- df_salience %>% 
  group_by(attribute, level) %>% 
  summarize(average = mean(selected)) %>%
  ungroup() %>% 
  mutate(deviation = abs(average - 0.5)) %>% 
  group_by(attribute) %>% 
  summarize(salience_mean = mean(deviation))

## Visualize the estimates
salience$attribute <- 
  factor(salience$attribute,
         levels = c("Leader\nSelection",
                    "Civil\nLiberties",
                    "Institutional\nChecks",
                    "National\nEconomy",
                    "Respondent\nWealth",
                    "Public\nSafety",
                    "Corruption\nin Politics",
                    "Health\nCare",
                    "Minority\nTreatment",
                    "Respondent\nIdentity"))
p_salience <- ggplot(data = salience) +
  geom_point(aes(x = fct_rev(attribute), y = salience_mean),
             size = 2, shape = 8) + 
  coord_flip(ylim = c(0, 0.08)) +
  labs(y = "Attribute salience", x = NULL) +
  theme_bw(base_size = 10, base_family = "Times") %+replace%
  theme(
    axis.text.x = element_text(size = 9, colour = "black", hjust = .5, vjust = 1),
    axis.text.y = element_text(size = 9, colour = "black", face = "bold", hjust = 1, vjust = .5),
    axis.ticks = element_line(colour = "grey50"),
    strip.text.x = element_text(size = 9),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()
  )

## Combine the plot with the aggregate AMCE plot
p_amce <- p_amce + ggtitle("Estimated AMCEs") +
  theme(plot.title = element_text(hjust = 0.5, size = 11, face = "italic"))
p_salience <- p_salience + ggtitle("Attribute Salience") + 
  theme(plot.title = element_text(hjust = 0.5, size = 11, face = "italic"))
plot_grid(p_amce, p_salience, labels = "AUTO", nrow = 1, 
          label_fontfamily = "Times", label_size = 11,
          rel_widths = c(2, 1))
ggsave("Figure S7.pdf", width = 8, height = 6)
